random number generation and thread safety

[vgoudar]

Im working on creating a point process and a membrane mechanism with NMODL and both of these require the use of random numbers. My membrane mechanism uses scop_random() and my point process uses set_seed() (similar to the example pregen.mod) and poisrand(). The NMODL compiler complains that all these are not thread safe, but surely there is a way to get them to be. A similar warning is shown for the example pregen.mod. I believe the correct fix is to add the THREADSAFE keyword. Am I right? Or is there a thread safe version of this functionality available?

[hines]

scop_random is not threadsafe because it uses global state which changes on each pick of a random number.

I much prefer the use of the Random class with MCellRan4 as used in nrn/src/nrnoc/netstim.mod. Point processes that each use their own distinct Random instance are not only threadsafe but also make it much easier to get reproducibility in parallel cluster computing.

With regard to density membrane mechanisms it will be sufficient to use a distinct Random instance per cell (different cells can be in different threads but all the compartments of any given cell are in the same thread (assuming no splitting)).

By the way, it is easier these days to get statistically independent random streams from the Random class because the low index for MCellRan4 is settable for each instance.
http://www.neuron.yale.edu/neuron/stati ... #MCellRan4

[vgoudar]

Im a bit new no neuron, so Im a little slow to it ... the Random and MCellRan4 classes seem applicable to hoc ... how would I then use them in nmodl? The netstim code uses nrn_random_pick and Im not sure how and nrn_random_arg and Im not sure how these are related to MCellRan4. Also, in the netstim code, what is the relevance of donotuse and _p_donotuse? Thanks for your help!

[hines]

netstim.mod has a noiseFromRandom() procedure which is called from user hoc code and takes an objref referencing a Random instance (which should be using the negexp distribution and the MCellRan4 generator). Think of that implementation as a pattern for you to use in linking your hoc Random instance and an instance of a point process.

nrn_random_pick(_p_donotuse) picks a random number from the hoc Random instance.

nrn_random_arg(1) gets the first argument (which must be a hoc object reference to a Random instance).

Code: Select all

void** pv = (void**)(&_p_donotuse);

is bit of magic that allows a pointer to a double to in fact keep a pointer to the Object* (the reference to the Random instance).

[vgoudar]

Thanks for your help! Here is my template:

Code: Select all

begintemplate ExSoma		// create a new template object...
public soma

objref nclistEx, nclistInh, noiseRandObj 

create soma
proc init() { 
    nclistEx = new List() // Instantiate list of netcons
    nclistInh = new List() // Instantiate list of netcons

    noiseRandObj = new Random() //Provides NOISE with random stream 
    noiseRandObj.uniform(0,1)
    
    soma {
	Ra = 35		// Axial Resistivity
	nseg = 1
	diam = 10	// diameter in um: L and diam are used to calculate area
	L = 10		// length in um: relevant because g_pas is in mho/cm2 
	cm = 3		// capacitance
	
	insert ExIAF		// includes
	ThrConst_ExIAF = -40		// -70 
	ePAS_ExIAF = -60	// -70
        gPAS_ExIAF = 0.0001	// time constant is 30 ms because cm = 3, 0.0001
        eOFF_ExIAF = -60	//e_pas
	gAHPbar_ExIAF = 0.00007
	//tauAHP_ExIAF  = 0.00000001
	insert NOISE
        setRandObjRef_NOISE(noiseRandObj)
    }
    //print "<< ExSoma >>"
}
endtemplate ExSoma

The NOISE membrane mechanism looks like this:

Code: Select all

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.uniform(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 setRandObjRef() {
VERBATIM
	void** pv4 = (void**)(&_p_randObjPtr);
	if (ifarg(1)) {
		*pv4 = nrn_random_arg(1);
	}else{
		*pv4 = (void*)0;
	}
ENDVERBATIM
}

When I then instantiate ExSoma from a file NETWORK.oc, I get the following error:

Code: Select all

/Applications/NEURON-7.1/nrn/umac/bin/nrniv.app/Contents/MacOS/nrniv: Bus error See $NEURONHOME/lib/help/oc.help
 in NETWORK.oc near line 81
 }
  ^
        ExSoma[0].setRandObjRef_NOISE(Random[2])
      ExSoma[0].init()
    xopen("NETWORK.oc")

What am I doing wrong?

[hines]

Calling functions in density mechansms which use range variables require that one
specify which instance of the mechanism is being referred to (ie. not only the section
but also the segment). So , assuming all your sections only have nseg=1,
add a line to the fragment

Code: Select all

   insert NOISE
        setRandObjRef_NOISE(noiseRandObj)

so it looks like

Code: Select all

   insert NOISE
        setdata_NOISE(0.5)
        setRandObjRef_NOISE(noiseRandObj)

Note this strategem is not necessary for POINT_PROCESS objects, only for SUFFIX mechanisms.