Making Random Processes Threadsafe

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

Moderator: hines

Post Reply
agmccrei
Posts: 24
Joined: Thu Sep 10, 2015 1:34 pm

Making Random Processes Threadsafe

Post by agmccrei »

Hello,

I have a question about using mod files that utilize random number generators in parallelized network simulations. Specifically I use Gfluct.mod (https://senselab.med.yale.edu/modeldb/S ... mod#tabs-2 - from what I've read, this file also appears to have been corrected to not have negative conductance values), which uses normrand(0,1). According to viewtopic.php?t=3008, this is not threadsafe to use in parallel simulations, but I'm not sure what the effect would be on a parallelized network simulation. Moreover, when I assign each Gfluct process on each neuron with it's own unique seed value (i.e. new_seed()), each neuron appears to have unique voltage and spiking - so not sure if it's even worth it to adapt Gfluct.mod to accept Random streams (i.e. as in https://neuron.yale.edu/phpBB/viewtopic ... 986#p12288).

In any case, I adapted Gfluct.mod to accept Random streams, as in ingauss.mod from the above link posted by Ted, and this seems to work when I run it and generates similar results:

Code: Select all

INDEPENDENT {t FROM 0 TO 1 WITH 1 (ms)}

NEURON {
	POINT_PROCESS Gfluct2
	RANGE g_e, g_i, E_e, E_i, g_e0, g_i0, g_e1, g_i1
	RANGE std_e, std_i, tau_e, tau_i, D_e, D_i
	RANGE new_seed
	NONSPECIFIC_CURRENT i
	THREADSAFE
    POINTER donotuse
}

UNITS {
	(nA) = (nanoamp) 
	(mV) = (millivolt)
	(umho) = (micromho)
}

PARAMETER {
	dt		(ms)

	E_e	= 0 	(mV)	: reversal potential of excitatory conductance
	E_i	= -75 	(mV)	: reversal potential of inhibitory conductance

	g_e0	= 0.0121 (umho)	: average excitatory conductance
	g_i0	= 0.0573 (umho)	: average inhibitory conductance

	std_e	= 0.0030 (umho)	: standard dev of excitatory conductance
	std_i	= 0.0066 (umho)	: standard dev of inhibitory conductance

	tau_e	= 2.728	(ms)	: time constant of excitatory conductance
	tau_i	= 10.49	(ms)	: time constant of inhibitory conductance
}

ASSIGNED {
	v	(mV)		: membrane voltage
	i 	(nA)		: fluctuating current
	g_e	(umho)		: total excitatory conductance
	g_i	(umho)		: total inhibitory conductance
	g_e1	(umho)		: fluctuating excitatory conductance
	g_i1	(umho)		: fluctuating inhibitory conductance
	D_e	(umho umho /ms) : excitatory diffusion coefficient
	D_i	(umho umho /ms) : inhibitory diffusion coefficient
	exp_e
	exp_i
	amp_e	(umho)
	amp_i	(umho)
    donotuse
}

INITIAL {
	g_e1 = 0
	g_i1 = 0
	if(tau_e != 0) {
		D_e = 2 * std_e * std_e / tau_e
		exp_e = exp(-dt/tau_e)
		amp_e = std_e * sqrt( (1-exp(-2*dt/tau_e)) )
	}
	if(tau_i != 0) {
		D_i = 2 * std_i * std_i / tau_i
		exp_i = exp(-dt/tau_i)
		amp_i = std_i * sqrt( (1-exp(-2*dt/tau_i)) )
	}
}

BREAKPOINT {
	SOLVE oup
	if(tau_e==0) {
	   g_e = std_e * grand()
	}
	if(tau_i==0) {
	   g_i = std_i * grand()
	}
	g_e = g_e0 + g_e1
	if(g_e < 0) { g_e = 0 }
	g_i = g_i0 + g_i1
	if(g_i < 0) { g_i = 0 }
	i = g_e * (v - E_e) + g_i * (v - E_i)
}


PROCEDURE oup() {		: use grand()
   if(tau_e!=0) {
	g_e1 =  exp_e * g_e1 + amp_e * grand()
   }
   if(tau_i!=0) {
	g_i1 =  exp_i * g_i1 + amp_i * grand()
   }
}


PROCEDURE new_seed(seed) {		: procedure to set the seed
	set_seed(seed)
	VERBATIM
	  printf("Setting random generator with seed = %g\n", _lseed);
	ENDVERBATIM
}

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

FUNCTION grand() {
VERBATIM
    if (_p_donotuse) {
        /*
         : Supports separate independent but reproducible streams for
         : each instance. However, the corresponding hoc Random
         : distribution MUST be set to Random.uniform(0,1)
         */
            _lgrand = nrn_random_pick(_p_donotuse);
    }else{
        /* only can be used in main thread */
        if (_nt != nrn_threads) {
hoc_execerror("multithread random in InUnif"," 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
        grand = normrand(0,1)
VERBATIM
    }
ENDVERBATIM
}

PROCEDURE noiseFromRandom() {
VERBATIM
 {
    void** pv = (void**)(&_p_donotuse);
    if (ifarg(1)) {
        *pv = nrn_random_arg(1);
    }else{
        *pv = (void*)0;
    }
 }
ENDVERBATIM
}

And the hoc code looks something like this (i.e. this is just a simplified version for illustrative purposes):

Code: Select all

Nmc = 150
N_gflucts_per_neuron = 10
for (neuron_number=0; neuron_number<Nmc; neuron_number+=1){
	for (gfluct_num=0; gfluct_num < N_gflucts_per_neuron; gfluct_num+=1){
		access treename[secNum]
		treename[secNum] sref = new SectionRef()
		sref {
			OUprocess.append(new Gfluct2(siteX))
			OUni = OUprocess.count()-1 // OU object index
			
			// Set OU parameters
			//
			// OUprocess.o[OUni].new_seed(neuron_number*N_gflucts_per_neuron+gfluct_num) // This appears to not be threadsafe
	  	 	rslist.append(new Random(neuron_number*N_gflucts_per_neuron+gfluct_num))
	  	 	rslist.o(OUni).normal(0,1)
	  	 	OUprocess.o(OUni).noiseFromRandom(rslist.o(OUni))
		}
	}
}
As I said, I'm not 100% sure if this is the right or wrong way of doing this type of thing in a parallelized network simulation, or if it's equivalent to just using normrand(0,1). Any clarity on the topic would be helpful.

Thanks,

Alex GM
Post Reply