proposal for a native random123 type in nmodl

[ohm314]

Random123 variables could be treated in NEURON as first class citizen. A modfile using such random variables could look like this:

NEURON	{ 
  ARTIFICIAL_CELL NetStim
  RANGE interval, number, start
  RANGE noise
  RANDOM123 r
}

PARAMETER {
	interval	= 10 (ms) <1e-9,1e9>: time between spikes (msec)
	number	= 10 <0,1e9>	: number of spikes (independent of noise)
	start		= 50 (ms)	: start of first spike
	noise		= 0 <0,1>	: amount of randomness (0.0 - 1.0)
}

ASSIGNED {
	event (ms)
	on
	ispike
	donotuse
}


INITIAL {

	setseq(r, 0, 0)

	on = 0 : off
	ispike = 0
	if (noise < 0) {
		noise = 0
	}
	if (noise > 1) {
		noise = 1
	}
	if (start >= 0 && number > 0) {
		on = 1
		: randomize the first spike so on average it occurs at
		: start + noise*interval
		event = start + invl(interval) - interval*(1. - noise)
		: but not earlier than 0
		if (event < 0) {
			event = 0
		}
		net_send(event, 3)
	}
}	

PROCEDURE init_sequence(t(ms)) {
	if (number > 0) {
		on = 1
		event = 0
		ispike = 0
	}
}

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*erand()
	}
}


FUNCTION erand() {
        erand = exprand(r)
}


PROCEDURE noiseFromRandom123(i, j, k) {
        rstream(r, i, j, k)
}


PROCEDURE next_invl() {
	if (number > 0) {
		event = invl(interval)
	}
	if (ispike >= number) {
		on = 0
	}
}

NET_RECEIVE (w) {
	if (flag == 0) { : external event
		if (w > 0 && on == 0) { : turn on spike sequence
			: but not if a netsend is on the queue
			init_sequence(t)
			: randomize the first spike so on average it occurs at
			: noise*interval (most likely interval is always 0)
			next_invl()
			event = event - interval*(1. - noise)
			net_send(event, 1)
		}else if (w < 0) { : turn off spiking definitively
			on = 0
		}
	}
	if (flag == 3) { : from INITIAL
		if (on == 1) { : but ignore if turned off by external event
			init_sequence(t)
			net_send(0, 1)
		}
	}
	if (flag == 1 && on == 1) {
		ispike = ispike + 1
		net_event(t)
		next_invl()
		if (on == 1) {
			net_send(event, 1)
		}
	}
}

[pramodk]

@halfflat : here is the relevant ticket to your question during CNS. Above is just draft from first discussion a year ago. We haven't started any implementation aspects yet but if you have experiences from your implementation and usage, it will be great to have your inputs.

[ramcdougal]

One question: what responsibility do we have to engage the broader community on the NMODL language description now that other tools (e.g. Arbor) use NMODL?

[pramodk]

In my previous comment I cc'd Sam Yates who is working with Arbor and knows NMODL/Random123 aspects very well. (@bcumming : tagging you as well in case Sam is away during summer vacation).

This item is on Friday's meeting agenda so hopefully we will discuss it there.

[ramcdougal]

Oh that's great. Thanks.

Another thing to consider, should it be possible to say GLOBAL RANDOM123 or RANGE RANDOM123 to indicate if it should be shared or per-instance? I guess the default would be per-instance? I think that's what happens with POINTER now?

[nrnhines]

I don't have a use case for GLOBAL in this context. However it does raise the question about default id1,id2, id3 stream indices and initialization. Naively, 0,0,0 will make all instances generate the same stream, and initialization when the INITIAL block should be at the choice of the user. In other words, perhaps there should be an error or warning if the user does not at least call setseq.

Also, setseq(r, ...) fits into standard NMODL function syntax, but it could be made to look more object oriented.

[pramodk]

Another question came to discussion yesterday : Do we need name this interface as RANDOM123? (i.e. 123 part could be removed). Random123 is just one implementation of RNG.

[halfflat]

Thanks for including me in the conversation!

We don't yet have a mechanism for specifying randomness in the NMODL that we parse in Arbor, though we will need one shortly for the case of modelling a stochastic ODE.

We have though of course been discussing this among ourselves, and we're likely to pursue a less explicit approach, where randomness is included by introducing terms describing stochastic processes. Random123 then is an implementation matter.

An aspirational goal is reproducibility across different distributions of the model across ranks and groups; we correspondingly aim to tie each stochastic source to a pseudo-random sequence that is determined entirely by its 'location' in the model definition, together with any run-time seed data given at simulation time. Random123 is naturally a really good fit for this sort of approach.

[ohm314]

I wanted to pick up this topic again.. In the meantime we also had a discussion in the dev meeting (here).

Two main questions I am wondering about:

• It makes sense to make this library agnostic (so remove the 123 from the type name). But the problem is that random123 has a specific API on how to setup seeds and streams requiring specific parameters (seed(r, i, j, k) and setseq(r, a, b)), other random number generators might have again other requirements. How can we hide the random123 specificities but still allow the user to initialize the streams correctly with different streams having different seeds? On the other hand if we ever support a serial RNG, there the notion of multiple streams wouldn't make sense and a simpler API with seed(r, s) would suffice.
• Should the API support random number distributions? For example in STL this is part of the random number specification. One initializes a generator, and using the generator one initializes a distribution from which one can draw. I think it would be useful to have this too. If we are going to add this language element, I think it would be good to think of this aspect too, rather than coming back to it later again.
• And in response to @nrnhines comment above, yes we can also use a more object oriented syntax (r.seed(i, j, k) and r.setseq(a, b)).

@nrnhines , @pramodk (and others) any thoughts?

[nrnhines]

I think a special meeting would be useful. Looking at slide 11 of the above dev meeting link. I'm indecisive about whether seed and setseq need to be explicit as opposed to r.start() which itself I can imagine the user may sometimes want (to reproduce an earlier sim) or not want (just continue where it last left off). At any rate it doesn't seem to be relevant in the mod file itself. On the other hand, specification of the distribution does seem relevant within the mod file. I'm also on the fence with regard to whether Random should be in the NEURON block or a block of its own. Not declared in the NEURON block if the random variable is not really managed by the interpreter. @halfflat mentions an intrigueing aspiration: Declaration only of the name and distribution of a random variable in the mod file and everything else is behind the scenes.

[ohm314]

After having met with Michael and other members of the dev team, here is how we think this new language feature should look like:

Declaration

• the random type should not be tied to an implementation (random123 or other) but rather simply express that the declared variable is sampled from a distribution function.
• Since it is a variable declaration, one logical place to put it is in the NEURON block.
• Random variables are all implicit RANGE variables, which means that for each instance of the mechanism one instance of the random variable is maintained.
• Syntactically this declaration would look as follows:

RANDOM DISTRIBUTION(param_1, ..., param_n) VARNAME_1, ..., VARNAME_M

where:

• RANDOM is a new reserved word
• DISTRIBUTION would be the actual name of the distribution function
  • Initially UNIFORM, EXP could be supported
• param_1, ..., param_n are the distribution function parameters, for example
  • UNIFORM(a, b) would define a uniformly distributed random number over the range (a,b).
  • EXP(lambda) would define an exponentially distributed random number with parameter lamba.

Example variable definition:

RANDOM UNIFORM(0, 1) ur1, ur2, ur3

Three uniformly i.i.d random variables over the range (0,1)

RANDOM EXP(2.0) er1

An exponetially distributed random number with parameter lambda=2.0.

Use

The random variables can be used in any block of the NMODL file (e.g. BREAKPOINT, PROCEDURE, ...). Using such a variable in an expression means that we sample from the distribution set at the declaration. Furthermore, sampling a random variable twice within the same time step should result in the same value. This is especially important in the context of BREAKPOINT where the same code block is called several times in the final generated code, and resampling the random variable at each execution of the block would result in inconsistent results.

Questions

• Should the default be that the seed is fixed to a constant value to allow exact reproducibility or should the default be that the seed is different at each run? The latter would mean that the random number sequence is not the same at two different runs but the random number distribution is of course still reproducible.
  • both behaviors seem useful, so enabling one or the other could be done with an extra keyword at variable declaration or through interpreter API.

I would really appreciate more feedback (if any) so we can nail down some more details and try out a prototype implementation.

[ramcdougal]

Sounds good.

Some questions:

• Are the values ur1, ur2, etc chosen from the distribution but themselves regular doubles? Alternatively, are they different every time they're accessed?
• If they're regular doubles and they are declared inside the NEURON block, is there a way to use them to generate some sort of irregular behavior at different times? (i.e. can we resample continuously throughout the simulation?)
• should it be possible to set the seeds separately for the different variables? I guess here I'm assuming the values are different each time they're accessed?

[nrnhines]

Since it is a variable declaration, one logical place to put it is in the NEURON block.

That is correct, but it is also the case that the NEURON block generally means to me what the model wishes to expose to NEURON. Ie. should the random variable be declared as GLOBAL or RANGE as for PARAMETER and ASSIGNED. In that case
the random_decl: RANDOM distribution params namelist production would appear as a top-level statement.
But I would accept it as a NEURON block declaration with the stipulation that it is always RANGE.

Are the parameters constants, or can they be the names of parameters. In practice, the negexp distribution always had the parameter 1.0 and the user multiplied the picked value by lambda (which had units of (ms)). Similarly it was generally the case that for the normal distribution, the mean was 0.0 and the variance was 1.0 and the picked value would be x*var + mean, again with consistent units. Not all distributions have this normalization property, e.g. poisson but I hate to drag units into the discussion if it would be possible to say all parameters and pick results are dimensionless.

I would vote for default "seed is fixed to a constant value to allow exact reproducibility" but also my opinion is that it is not the business of the model description to specify that kind of higher level user specifiable property I.e. in what way is the random variable not random. Sometimes runs are looking for reproducibility for parallel debugging and sometimes runs are collecting statistically independent results. I guess I would change my vote to say that it is a higher level simulator property and subject to implementation feasibility. Nothing about it can be specified in nmodl and nothing about it is implicit in the nmodl translation.
From nmodl's point of view x.pick is a complete black box except with respect to the distribution chosen.

In my opinion a sample implementation can begin and end with Random123 and it suffices to allow distributions of uniform, negexp, and normal. I'd go so far as to say the user can leave out the () and expect the normalized distribution.

sampling a random variable twice within the same time step should result in the same value.

This is nice if it is implementable. With respect to cvode, the only safe place to do it is in a BEFORE STEP block. It may be satisfactory to follow the policy of Random.play . See https://neuronsimulator.github.io/nrn/py_doc/programming/math/random.html#Random.play which, in a kind of backhanded
incomplete way is an implementation of this proposal except that in the mod file the random variable is neither declared nor is its distribution specified. With respect to random variables in NET_RECEIVE block the user would expect a new value on every evaluation of the variable (in concert with its event driven nature). If a random variable appears in a function, it gets a little tricky without knowing the caller that determines whether evaluation picks a new value.

[ramcdougal]

Backing up a step and speaking for nobody other than myself...

My original conception of this was as a way of eliminating a class of VERBATIM block usage by standardizing random specification. I still think this, by itself, would be progress.

I love the idea of random variables that know to resample themselves at the beginning of each timestep and then act as constants throughout. That said, I cannot of any analogous case of specifications that have values that change randomly without an explicit assignment, and I could see this as a surprise like you'd get from a global variable that isn't set anywhere you can see.

In XPPAUT, an equation could directly invoke e.g. normal(0,1) as in this example:

ie'=(-ie+a1*normal(0,1))*w1

Since order of equations mostly doesn't matter in XPPAUT, the same number could presumably be reused within a timestep by explicitly assigning it to a variable. We could do the same.

The upshot of a function-based approach is that it doesn't require any changes to the parser to allow picking from a (presumably) location-aware random stream. Ideally, each instance pulls from its own stream so that we don't have to worry about what is evaluated in what order (thus avoiding any sort of procedural programming type side effects).

Another option that I like less, but just to have it on the table, would be to introduce a bit of object-orientedness into NMODL and allow a Random class that's essentially the same thing as NEURON's Random class. Weak objection: it's possibly tied too closely to NEURON and it's not clear what this would mean in e.g. Arbor. Stronger objection: I'm not sure if this makes NMODL farther away from a declarative format but I don't think it helps.

Finally, I think @nrnhines 's point about NET_RECEIVE calling functions vs other function calls is something that a solution would have to address.

[ohm314]

Thanks a lot for the feedback!

That is correct, but it is also the case that the NEURON block generally means to me what the model wishes to expose to NEURON. Ie. should the random variable be declared as GLOBAL or RANGE as for PARAMETER and ASSIGNED. In that case the random_decl: RANDOM distribution params namelist production would appear as a top-level statement.
But I would accept it as a NEURON block declaration with the stipulation that it is always RANGE.

Yes, these variables would be implicitly always RANGE, and now that you mention it, actually it does make sense that they are in the NEURON block, because they should be exposed the NEURON as we want to be able to do advanced / user-defined manipulations via the interpreter interface. The user doesn't have to do anything with those variables from the NEURON side, but they can if they wish to (e.g. to set the seed)

Are the parameters constants, or can they be the names of parameters.

I think it should be possible to use constants or parameters. To me this makes the whole random distribution concept much more complete and clean if I can set the distribution with it's parameters at declaration and don't ahve to worry about applying the transformation myself whenever I sample. I "think" it shouldn't be much harder to provide this as part of the specification and it would clean up the NMODL code some more by making such transformations unnecessary.

I would vote for default "seed is fixed to a constant value to allow exact reproducibility" but also my opinion is that it is not the business of the model description to specify that kind of higher level user specifiable propert

Absolutely, but we still need to pick one behavior that happens by default. The behavior is not changed in the NMODL code because as you said that is business of the simulator and not the model description. But if in the simulator we don't set anything either, then implicitly we assume that the random number seed is fixed (or not...).

In my opinion a sample implementation can begin and end with Random123 and it suffices to allow distributions of uniform, negexp, and normal. I'd go so far as to say the user can leave out the () and expect the normalized distribution.

OK, cool, then we begin with uniform, negexp and normal as distributions. I'm wondering if we should make the () optional. The thing is, when this is optional, then the user will have to look up the defaults, arguably it's quite obvious here, but maybe safer to simply ask the user to set them and then it's written down right in the model.

sampling a random variable twice within the same time step should result in the same value.

Regarding this: After reading @nrnhines and @ramcdougal 's comments, I think that as a concept it would be nice to have this, but in practice it will be complicated to make it work correctly. Additionally as @ramcdougal points out, this may surprise users quite a bit, since it is not something typically seen. I would say we better keep it simple. Whenever a random variable is accessed, a new value is sampled from the distribution. If we would like to reuse a value to do multiple computations (e.g. in multiple CVODE steps) then we must first do a x = r and continue with x.

The idea of using a function syntax instead is interesting too, it would certainly simplify things for us, but at the same time, I find it will be quite wordy, and I sort of like the neat/clean code that random variables would introduce. Would you agree that with above simplification (using the variable always means we are resampling) we get a simple syntax without the semantic ambiguity and complications? It also shouldn't be a problem for the parser, or am I missing something?

Thanks again @nrnhines and @ramcdougal for your input!

[nrnhines]

because they should be exposed the NEURON as we want to be able to do advanced / user-defined manipulations via the interpreter interface.

That reminds me of an issue I didn't think very consistently about. What does access to the random variable even mean in the interpreter? We do not want it to mean the pick of a new random value. I was thinking in terms of properties such as initialization of the stream (and stream identifiers). And setting distribution parameters for the random variable. That seems very much like an object. But that seems analogous to the
hoc Random object and close to @ramcdougal original interpretation of this project devoted primarily to eliminating the VERBATIM block by introducing an NMODL aware Random object declaration. I still think the location default to set stream id is really useful as it offloads a great deal of cognitive effort from the user. And from @ramcdougal interpretation it would mean the automatic creation of a hoc Random instance with properties such as h.Random123(id1, id2, id3) and h.Random.negexp(1). I guess I was imagining a completely new python interface object that would know about the hidden internal instance of nrnran123_State which would be, in fact, the distinct structure associated with each the NMODL RANDOM``` statement variable.

There is the edge case of NetStim without a gid in a parallel sim where the NetStim must remain on the same process as its targets (and for CoreNEURON, in the same thread). That complicates the automatic location choice of stream id (often determined by the (single) target the NetStim is connected to (that does have a gid)).

[ohm314]

We do not want it to mean the pick of a new random value. I was thinking in terms of properties such as initialization of the stream (and stream identifiers). And setting distribution parameters for the random variable. That seems very much like an object.

Exactly my thinking. The interpreter access the object that is controlling the random variables in the model

And from @ramcdougal interpretation it would mean the automatic creation of a hoc Random instance with properties such as h.Random123(id1, id2, id3) and h.Random.negexp(1)

Not sure I understand this exactly.. do you mean that objects can be created from the interpreter without being declared in the model too? Where would the "live"? Not even sure if my question makes sense.

I guess I was imagining a completely new python interface object that would know about the hidden internal instance of nrnran123_State which would be, in fact, the distinct structure associated with each the NMODL RANDOM``` statement variable.

To me one RANDOM variable, as it is implicitly RANGE has as many instances as there are mechanism instances. The same variable is accessible through the interpreter, but the individual instances are not accessible. What is exposed is the object that holds the generator and the distribution from which the instances sample. Does that make sense?

[nrnhines]

individual instances are not accessible

As far as that can be carried out, I prefer it. But there are instance specific properties such as the stream id triple and (of less importance) the current sequence value. And the distribution parameters? The latter two will almost 100% of the time be the same for all instances. ie. By default the current sequence value set to 0 on initialization for all instances (of a specific random variable type). Typical access to mod file variables are on a per instance basis (RANGE as opposed to GLOBAL). This does not
necessarily prejudice against a fancy python object that holds the entire per process collection and at least allows type wide property setting.

And now I can't remember. Were we going to do the same thing for hoc Vector? ie. language support to avoid VERBATIM

do you mean that objects can be created from the interpreter without being declared in the model too?

Interpreter objects can be created from mod file statements (in VERBATIM blocks). The relevant interface for hoc Vector is

extern "C" {
extern Vect* vector_new(int, Object*); // use this if possible
extern Vect* vector_new0();
extern Vect* vector_new1(int);
extern Vect* vector_new2(Vect*);
extern void vector_delete(Vect*);
extern int vector_buffer_size(Vect*);
extern int vector_capacity(Vect*);
extern void vector_resize(Vect*, int);
extern Object** vector_temp_objvar(Vect*);
extern double* vector_vec(Vect*);
extern Object** vector_pobj(Vect*);
extern Vect* vector_arg(int);
extern int is_vector_arg(int);
extern char* vector_get_label(Vect*);
extern void vector_set_label(Vect*, char*);
}

Where in mod VERBATIM Vect* and Object* are replaced by void*. Hoc Random has a few of these also but repick, ids, and set_sequence are the only ones I recall and I don't think the constructor is available.

[ramcdougal]

If I change nseg on the apical dendrites, would that change the random streams used on the basals? (It'd be best if not.)

[nrnhines]

If I change nseg on the apical dendrites, would that change the random streams used on the basals?

I think the algorithm for converting (cell, section, x, mech, name) -> (id1, id2, id3) should preserve already existing stream id's as much as possible.

[bcumming]

I have an update from Arbor, now that we implemented support for stochastic differential equations.

To start, I will refer to some pages in the Arbor docs.

General documentation about support for SDEs in Arbor:
https://docs.arbor-sim.org/en/latest/dev/sde.html#stochastic-differential-equations

An illustrated example of its use in a single cell model:
https://docs.arbor-sim.org/en/latest/tutorial/calcium_stdp_curve.html

Documentation of the NMODL extension that supports this feature:
https://docs.arbor-sim.org/en/latest/fileformat/nmodl.html#stochastic-processes

This has been implemented such that a model will generate the same results between runs with the same seed, and it will produce the same results if the number of threads/MPI ranks is varied between runs, and if GPUs or CPU back ends are used.

The feature was implemented by Fabian, @boeschf.

[ohm314]

Just to keep the different pieces of information linked in one place. Here's a gist I had put together a while back with some example generated code:
https://gist.github.com/ohm314/1d47a6cfecc6f71d238161358aaf59d1

[ohm314]

Update on syntax

With @pramodk I've looked again over our example code and we found a few improvements (and worked out some of the questions that we had still remaining). I'd like to share with you now the updated proposal of how such variables could look like and how the generated code (with NMODL), along with the python interface, could look like. For simplicity I've abbreviated the netstim.mod example to only contain the parts that really matter to us here.

NEURON  {
  ARTIFICIAL_CELL NetStim
  RANGE interval, number, start
  RANGE noise
  THREADSAFE : only true if every instance has its own distinct Random
  : !!! this is how we declare random variables, not much has changed here. The declaration parsing is already implemented.
  RANDOM NEGEXP(1.0) r
}


INITIAL {

  : !!! there is no way to express this in a way that is not specific to Random123
  : we should think about whether we want to have Random123 APIs inside nmodl, or whether such an
  : "advanced" use-case is so rarely used that it can remain a C call in VERBATIM.
  setseq(r, 0, 0)

}

FUNCTION invl(mean (ms)) (ms) {
  if (mean <= 0.) {
    mean = .01 (ms) : I would worry if it were 0.
  }
  if (noise == 0) {
    invl = mean
  }else{
    : !!! we simply use the variable r. Whenever it is found in an expression, the distribution is sampled
    : NB. assigning to r is an error
    : If the same sampled random value is to be reused, first assign it to a local variable
    : LOCAL x
    : x = sample(r)
    invl = (1. - noise)*mean + noise*mean*sample(r)
  }
}

The generated C++ code looks like that:

#include <math.h>



#include <coreneuron/utils/randoms/nrnran123.h>
// include all the rest of headers


namespace coreneuron {

    /** all mechanism instance variables */
    struct NetStim_Instance  {
        const double* __restrict__ interval;
        const double* __restrict__ number;
        const double* __restrict__ start;
        double* __restrict__ noise;
        double* __restrict__ event;
        double* __restrict__ on;
        double* __restrict__ ispike;
        double* __restrict__ v_unused;
        double* __restrict__ tsave;
        const double* __restrict__ node_area;
        // !!! generated for each RANDOM variable to hold the Random123 counter struct
        void** __restrict__ rndcntr_r;
        void** __restrict__ point_process;
        void** __restrict__ tqitem;
    };


    /** initialize mechanism instance variables */
    static inline void setup_instance(NrnThread* nt, Memb_list* ml)  {
        NetStim_Instance* inst = (NetStim_Instance*) mem_alloc(1, sizeof(NetStim_Instance));
        int pnodecount = ml->_nodecount_padded;
        Datum* indexes = ml->pdata;
        inst->interval = ml->data+0*pnodecount;
        inst->number = ml->data+1*pnodecount;
        inst->start = ml->data+2*pnodecount;
        inst->noise = ml->data+3*pnodecount;
        inst->event = ml->data+4*pnodecount;
        inst->on = ml->data+5*pnodecount;
        inst->ispike = ml->data+6*pnodecount;
        inst->v_unused = ml->data+7*pnodecount;
        inst->tsave = ml->data+8*pnodecount;
        inst->node_area = nt->_data;
        // !!!
        // generate for each random variable a hidden pointer to its random123 counter struct
        inst->rndcntr_r = nt->_vdata;
        // !!!
        inst->point_process = nt->_vdata;
        inst->tqitem = nt->_vdata;
        ml->instance = (void*) inst;
    }



    inline double invl_NetStim(int id, int pnodecount, NetStim_Instance* inst, double* data, const Datum* indexes, ThreadDatum* thread, NrnThread* nt, double v, double mean);

    // !!! auto-generated hook to call from hoc/py (cell.init_r(gid, id1, id2)
    // takes in hoc/py 2 or 3 parameters, this is random123 specific but we hide that fact
    inline int init_r_NetStim(int id, int pnodecount, NetStim_Instance* inst, double* data, const Datum* indexes, ThreadDatum* thread, NrnThread* nt, double v) {
        int ret_init_r = 0;
        nrnran123_State** pv = (nrnran123_State**)(&inst->rndcntr_r[indexes[2*pnodecount + id]]);

        if (*pv) {
            nrnran123_deletestream(*pv);
            *pv = nullptr;
        }
        if (ifarg(3)) {
            *pv = nrnran123_newstream3((uint32_t)*getarg(1), (uint32_t)*getarg(2), (uint32_t)*getarg(3));
        }else if (ifarg(2)) {
            *pv = nrnran123_newstream((uint32_t)*getarg(1), (uint32_t)*getarg(2));
        }

        return ret_init_r;
    }
    // !!!


    inline double invl_NetStim(int id, int pnodecount, NetStim_Instance* inst, double* data, const Datum* indexes, ThreadDatum* thread, NrnThread* nt, double v, double mean) {
        double ret_invl = 0.0;
        if (mean <= 0.0) {
            mean = .01;
        }
        if (inst->noise[id] == 0.0) {
            ret_invl = mean;
        } else {
            // !!!
            // generate the appropriate call to the random number function possibly
            // with the necessary transformation into the requested distribution
            // Alternatively, add a little helper function sample(r), which contain the actual call
            // to nrnran123.
            // !!!
            ret_invl = (1.0 - inst->noise[id]) * mean + inst->noise[id] * mean * nrnran123_negexp((nrnran123_State*)inst->rndcntr_r[indexes[2*pnodecount + id]])r;
        }
        return ret_invl;
    }



    /** initialize channel */
    void nrn_init_NetStim(NrnThread* nt, Memb_list* ml, int type) {
        int nodecount = ml->nodecount;
        int pnodecount = ml->_nodecount_padded;
        const int* __restrict__ node_index = ml->nodeindices;
        double* __restrict__ data = ml->data;
        const double* __restrict__ voltage = nt->_actual_v;
        Datum* __restrict__ indexes = ml->pdata;
        ThreadDatum* __restrict__ thread = ml->_thread;

        setup_instance(nt, ml);
        NetStim_Instance* __restrict__ inst = (NetStim_Instance*) ml->instance;

        if (_nrn_skip_initmodel == 0) {
            int start = 0;
            int end = nodecount;
            #pragma ivdep
            for (int id = start; id < end; id++) {
                inst->tsave[id] = -1e20;
                double v = 0.0;
                // !!! the NMODL Random123 API could be translated into this (or this comes from VERBATIM)
                nrnran123_setseq((nrnran123_State*)inst->rndcntr_r[indexes[2*pnodecount + id]], 0, 0);
                inst->on[id] = 0.0;
                inst->ispike[id] = 0.0;
           }
        }
    }
}

Finally, here's how RANDOM variables are initialized (ie. seeded):

from neuron import h
pc = h.ParallelContext()
h.dt = 1.0/32


def test_bbcore():
    from neuron import coreneuron
    coreneuron.enable = True

    ncell = 10
    cells = []
    gids = range(pc.id(), ncell, pc.nhost()) # round robin
    for gid in gids:
        pc.set_gid2node(gid, pc.id())
        cell = h.NetStim()
        # !!! new way of initializing the seed for any RANDOM variable
        cell.init_rng(cell.r, gid, 2, 3)
        cells.append(cell)

    pc.spike_record(-1, tvec, idvec)

    h.CVode().cache_efficient(1)
    pc.set_maxstep(10)
    h.finitialize(-65)
    pc.psolve(10)

if __name__ == "__main__":
    test_bbcore()

@nrnhines , @ramcdougal (and others) : What do you think?

[pramodk]

Just to update here: @ohm314 - me & @nrnhines went through the example. There was comment about syntax: sample(r) vs r.sample(). But that's minor and we can check that bit later. So all good for implementing this.

[ohm314]

From the discussion with @nrnhines : setseq(r) could also be called from the interpreter.
IN the interpeter r is an object that has public functions to

• setseq (setting the sequence id)
• init (setting up a new stream)
• sample