Noisy current using normrand?

General issues of interest both for network and
individual cell parallelization.

Moderator: hines

Post Reply
ChrisR

Noisy current using normrand?

Post by ChrisR »

Hello everybody!

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()        
Ifluct1.mod:

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
}
ted
Site Admin
Posts: 6289
Joined: Wed May 18, 2005 4:50 pm
Location: Yale University School of Medicine
Contact:

Re: Noisy current using normrand?

Post by ted »

ChrisR wrote: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.
Probably because every mechanism instance that needs randomness is drawing its numbers from the same pseudorandom sequence generator. If that generator produces the sequence
x0, x1, x2 . . .
and you have only one variable a0 that needs to show random fluctuations, that variable will have the same sequence of values i.e.
x0, x1, x2 . . .
Now add a second variable a1 that also needs to fluctuate. It's going to eat up some of the xi values that would otherwise have gone to a0. Add a third variable that must fluctuate, and you can see that it must steal some of the values that would have gone to a0 and a1.

You need to ensure that each variable that needs to fluctuate is driven by its own pseudorandom sequence generator. Furthermore, you will want each sequence generator to produce a stream of values that is statistically independent of the streams produced by all other sequence generators. You definitely do NOT want each generator to have the same seed. It would be good to read the discussion in
Randomness in NEURON models
which you will find among the links at
http://www.neuron.yale.edu/neuron/docs
The particular example discussed there involves spike streams in which spikes occur at pseudorandom times, but the solution is quite similar: you'll need a separate instance of the RandomStream class for every variable that is to show random fluctuations in the course of a simulation.
ChrisR

Re: Noisy current using normrand?

Post by ChrisR »

Thanks for your reply ted.

First of all using the same seed is NOT my goal. It's just a test to see weather I am able to control the random sequence generator.
And it seems that its not working.

So what actually is the problem here? Inferring from your reply I suspect that the nmod function normrand is picking from the same generator, and does not, as I would have expected, generate a new one for every instance of the point process.
Thus all the models using a derivative of the Gfluct.mod implementation just would not work on a parallel architecture?

Thus the only way to go is to work is to set up an external generator for every instance of the point process and pass the variables to the point process?
But how would I do that? Are there codes out there where this has been already accomplished?

Thanks
Christian
ChrisR

Re: Noisy current using normrand?

Post by ChrisR »

Ok, so, after some "scavenger hunting" I think I've figured it out! It's giving my reproduceable noise streams independent of the number of hosts.

I've been using information from here: http://www.neuron.yale.edu/neuron/node/60
And this post here: http://www.neuron.yale.edu/phpBB/viewto ... =31&t=2136

However a little comment from you guys on whether I did something reasonable here would be great!

Especially on the following:
Since I don't know how long my streams will be I am setting each generator "noiseRandObj" with a different unique lowindex which is in fact the gid of the cell:
noiseRandObj.MCellRan4(1, gidlist[l])
Is this reasonable? To what should I set the highindex? I guess it needs something > 0 so it is not set by NEURON itself to something depending on the number of hosts.

Thanks
Christian


Ifluct_test.py: (execute with "mpiexec -n X python Ifluct_test2.py" with X=1-3)

Code: Select all

from pylab import *
from mpi4py import MPI
from neuron import h

N = 3
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 = []
noises = []
for l in range(size(gidlist)):  # for every cell in the gidlist 

    ### NEW ###
    noiseRandObj = h.Random()  # provides NOISE with random stream
    noiseRandObj.MCellRan4(1, gidlist[l])  # set lowindex to gid, set highindex to what?   
    noiseRandObj.normal(0,1)
    noises.append(noiseRandObj)        
    ###########
        
    fluct = h.Ifluct1(cells[l](0.5))
    fluct.m = 0      # [nA]
    fluct.s = 1      # [nA]
    fluct.tau = 1    # [ms]
    flucts.append(fluct)   # add to list

    ### NEW ###    
    fluct.setRandObj(noiseRandObj)
    ###########
    
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(10)

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()        

Ifluct1.mod:

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
    THREADSAFE : only true if every instance has its own distinct Random
    POINTER randObjPtr
}


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
    randObjPtr
}

INITIAL {
    x = m               : to reduce the transient, the state is set to its (expected) steady-state
}


BEFORE BREAKPOINT {
    noise = randGen()
}

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 }
}

VERBATIM
double nrn_random_pick(void* r);
void* nrn_random_arg(int argpos);
ENDVERBATIM

FUNCTION randGen() {
VERBATIM
   if (_p_randObjPtr) {
      /*
      :Supports separate independent but reproducible streams for
      : each instance. However, the corresponding hoc Random
      : distribution MUST be set to Random.normal(0,1)
      */
      _lrandGen = nrn_random_pick(_p_randObjPtr);
   }else{
      hoc_execerror("Random object ref not set correctly for randObjPtr"," only via hoc Random");
   }
ENDVERBATIM
}

PROCEDURE setRandObj() {
VERBATIM
   void** pv4 = (void**)(&_p_randObjPtr);
   if (ifarg(1)) {
      *pv4 = nrn_random_arg(1);
   }else{
      *pv4 = (void*)0;
   }
ENDVERBATIM
}
ted
Site Admin
Posts: 6289
Joined: Wed May 18, 2005 4:50 pm
Location: Yale University School of Medicine
Contact:

Re: Noisy current using normrand?

Post by ted »

Thanks for digging into this on your own--I'm at the SFN meeting and have had little time to devote to the Forum.
ChrisR wrote:Since I don't know how long my streams will be I am setting each generator "noiseRandObj" with a different unique lowindex which is in fact the gid of the cell:
noiseRandObj.MCellRan4(1, gidlist[l])
That should work.
To what should I set the highindex? I guess it needs something > 0 so it is not set by NEURON itself to something depending on the number of hosts.
Pick any nonzero integer you like.

And, of course, it is a good idea to test what you have done on a toy model--something small that allows you to run quick "sanity check" simulations" that generate just enough results to enable verification that things are working properly. This is where modular program design, with parameters that specify model size & complexity, is handy.
Post Reply