Random Number Generator in NetStim

[nonfcomm]

What random number generator does NetStim use? If it is not of cryptographic quality, what would be the easiest way to change the random number generator?

Thank you

[ted]

All NetStims get their random intervals from the same generator. The default is "a variant
of the Linear Congruential Generator (algorithm M) described in Knuth, Art of Computer
Programming, Vol. III in combination with a Fibonacci Additive Congruential Generator" (see http://www.neuron.yale.edu/neuron/stati ... m.html#ACG).
However, you can specify a better generator with use_mcell_ran4() -- see
http://www.neuron.yale.edu/neuron/stati ... mcell_ran4
This cryptographic quality generator is described here:
http://www.neuron.yale.edu/neuron/stati ... #MCellRan4

The next question you will probably ask is "how do I set the seed?"
With either generator, just call NetStim.seed(x). However, no matter how many NetStims
you have, every one of them uses the same generator so there is only one seed.
Consequently you cannot vary only one stream in a population of NetStims.

"What if that's exactly what I need to do?"
To quote an email from Michael Hines dated 10/25/2004, "copy and modify the
nrn/src/nrnoc/netstim.mod file so that it made use of mcellran4 directly
with separate lowindex values."

[mflynn]

Please forgive me, but I've tried to change netstim.mod and any changes give me an error, such as just changing the name of the random function:

Code: Select all

FUNCTION invl(mean (ms)) (ms) {
	if (mean <= 0.) {
		mean = .01 (ms) : I would worry if it were 0.
	}
	if (noise == 0) {
		invl = mean
	}else{
		invl = (1. - noise)*mean + noise*mean*nwrand()
	}
}
VERBATIM
double nrn_random_pick(void* r);
void* nrn_random_arg(int argpos);
ENDVERBATIM

FUNCTION nwrand() {
VERBATIM
	if (_p_donotuse) {
		/*
		:Supports separate independent but reproducible streams for
		: each instance. However, the corresponding hoc Random
		: distribution MUST be set to Random.negexp(1)
		*/
		_lerand = nrn_random_pick(_p_donotuse);
	}else{
ENDVERBATIM
		: the old standby. Cannot use if reproducible parallel sim
		: independent of nhost or which host this instance is on
		: is desired, since each instance on this cpu draws from
		: the same stream
		nwrand = exprand(1)
VERBATIM
	}
ENDVERBATIM
}

Is erand the only thing you can name the function?
What I am trying to do is to change NetStim so that it uses a different random seed each time.

[ted]

What I am trying to do is to change NetStim so that it uses a different random seed each time.

Then don't tinker with the NetStim's NMODL source code. Use a hoc statement that changes the seed. Read about the Random class in the Programmer's Reference.

[mflynn]

I've read the programmer's reference, but I can't find anything about how to automatically change the seed (say by grabbing the least significant digits from the system clock). I guess I could just manually change it every time I run a simulation, but that seems kinda inelegant and cumbersome.

[mflynn]

ok, I've found a bit of code that seems to pick a new seed everytime you run a simulation:

Code: Select all

r = new Random(2232)
temp = r.uniform(0,1)
for i=0,i<iSeed {temp2 = r.repick}

what I would like to know is how to do this with the poisson process used by NetStim, without editing the NMODL code.

Mark

[ted]

what I would like to know is how to do this with the poisson process used by NetStim, without editing the NMODL code.

Why is this even an issue? Do it with hoc, so you have control over the value of the seed. Then you can reproduce any simulation run you want, by explicitly setting the seed. If you need to generate a sequence of simulations with different seeds, first create a text file that contains a sequence of different seeds. Read this text file into a Vector, and for each successive run, use the next element of the Vector as the seed.

[mflynn]

I guess my problem is I don't know how to do it with hoc. netstim.mod uses a random number generator to produce a poisson distributed train of spikes. How do I change the seed for that random number generator? I thought maybe since

Note that multiple instances of the Random class will produce different streams of random numbers only if their seeds are different.

just creating another instance of the Random class with a seed that I manually changed each simulation run would produce different random sequence for each simulation, but that doesn't happen. How do I change the seed for netstim?

Mark

[ted]

How do I change the seed for netstim?

objrefname.seed(seedvalue)

Example: a NetStim generates a sequence of events governed by negexp distribution. The event times are recorded to a Vector that is printed at the end of each run.

Code: Select all

load_file("nrngui.hoc") // load standard run system etc.

use_mcell_ran4(1) // make scop_random in NMODL use mcell_ran4 generator
SEED_0 = 178205489
SEED_1 = 73887633

objref r
r = new Random()

objref ns, tvec, nc, nil
ns = new NetStim()
ns.interval = 1
ns.start = 0
ns.noise = 1 // negexp
tvec = new Vector()
nc = new NetCon(ns, nil)
nc.record(tvec)
objref nc // source continues to record to tvec, even though nc is destroyed

tstop = 3 // long enough since ns.start==0, ns.interval==1 

proc test() {
  if ($1!=0) ns.seed($1)
  run()
  printf("run %d", i)
  tvec.printf("\t%g")
  printf("\n")
}

print "===================="
print "3 runs with one seed"
print "===================="
for i=0,2 test(SEED_0)
print "========================"
print "3 runs with another seed"
print "========================"
for i=0,2 test(SEED_1)
print "==========================="
print "3 runs with different seeds"
print "==========================="
for i=0,2 test(SEED_0 + i*SEED_1)
print "================================="
print "3 runs without resetting the seed"
print "================================="
for i=0,2 test(0)

This program generates the following output:

Code: Select all

====================
3 runs with one seed
====================
run 0   0.505711        1.52482 2.29216
run 1   0.505711        1.52482 2.29216
run 2   0.505711        1.52482 2.29216
========================
3 runs with another seed
========================
run 0   0.30512 1.95959 2.51157
run 1   0.30512 1.95959 2.51157
run 2   0.30512 1.95959 2.51157
===========================
3 runs with different seeds
===========================
run 0   0.505711        1.52482 2.29216
run 1   0.189026        0.923573        1.06228
run 2   0.0489857       1.1928  1.62423 2.29101
=================================
3 runs without resetting the seed
=================================
run 0   1.14853
run 1   0.27276 0.577038        0.747876        1.24377 2.45078
run 2   0.048435        1.13385 1.82803 2.3581

Two items are of note. First, see how use_mcell_ran4() was employed to specify that the MCellRan4 generator is used.

Second, a rationale for explicitly specifying the seed at each run:
Although calling run() without resetting the seed will indeed produce a different set of event times on each call, the event times are all part of the same sequence that was started with whatever seed was (explicitly or implicitly) used for the very first run. Consequently, the event times in each run after the first run will depend on the number of events generated in all previous runs. This is clearly bad for reproducibility--the results produced by any run should depend only on user-specified parameters that pertain to that run. They should not depend on "accidents" like how many runs were executed previously, how long any of those runs lasted, or how many events were generated in the course of previous runs. It is far better to explicitly specify a new seed before each run.

Another comment:

Looking at netstim.mod for NEURON 7.0 alpha, I see the following:

Code: Select all

NEURON	{ 
 . . .
  THREADSAFE : only true if every instance has its own distinct Random
  POINTER donotuse
}
 . . .
FUNCTION erand() {
VERBATIM
	if (_p_donotuse) {
		/*
		:Supports separate independent but reproducible streams for
		: each instance. However, the corresponding hoc Random
		: distribution MUST be set to Random.negexp(1)
		*/
		_lerand = nrn_random_pick(_p_donotuse);
	}else{
		/* only can be used in main thread */
		if (_nt != nrn_threads) {
hoc_execerror("multithread random in NetStim"," only via hoc Random");
		}
ENDVERBATIM
		: the old standby. Cannot use if reproducible parallel sim
		: independent of nhost or which host this instance is on
		: is desired, since each instance on this cpu draws from
		: the same stream
		erand = exprand(1)
VERBATIM
	}
ENDVERBATIM
}

which makes me wonder if it is now possible for each NetStim to sample from a different pseudorandom sequence.

[ted]

makes me wonder if it is now possible for each NetStim to sample from a different pseudorandom sequence.

But of course it is, and has been for some time. And I should have remembered this from our contribution to
Brette et al.
Simulation of networks of spiking neurons: a review of tools and strategies.
J. Comput. Neurosci. 23:349-398, 2007
(preprint available from http://www.neuron.yale.edu/neuron/bib/nrnpubs.html, code available from ModelDB via accession number 83319--http://senselab.med.yale.edu/modeldb/Sh ... odel=83319).
But somehow, outside of the context of distributed processing on parallel hardware, I completely forgot about it.

The strategy is to create, for each instance of the NetStim class, a corresponding instance of the Random class, making sure that the latter picks from negexp(1). Then use the NetStim's noiseFromRandom() PROCEDURE to link the two. For an example go to
http://senselab.med.yale.edu/modeldb/sh ... odel=83319
and drill down to the file
destexhe_benchmarks/NEURON/common/netstim.hoc

[mflynn]

Thanks! That did the trick.

Mark