On Tue, Mar 12, 2013 at 9:25 PM, Nathaniel Smith wrote:

On Mon, Mar 11, 2013 at 9:46 AM, Robert Kern wrote:

...

On Sun, Mar 10, 2013 at 6:12 PM, Siu Kwan Lam wrote:

...

My suggestion to overcome (1) and (2) is to allow the user to select between the two implementations (and possibly different algorithms in the future). If user does not provide a choice, we use the MT19937-32 by default.

numpy.random.set_state("MT19937_64", …) # choose the 64-bit implementation

Most likely, the different PRNGs should be different subclasses of RandomState. The module-level convenience API should probably be left alone. If you need to control the PRNG that you are using, you really need to be passing around a RandomState instance and not relying on reseeding the shared global instance.

+1

...

Aside: I really wish we hadn't exposed `set_state()` in the module API. It's an attractive nuisance.

And our own test suite is a serious offender in this regard, we have tests that fail if you run the test suite in a non-default order... https://github.com/numpy/numpy/issues/347

I wonder if we dare deprecate it? The whole idea of a global random state is just a bad one, like every other sort of global shared state. But it's one that's deeply baked into a lot of scientific programmers expectations about how APIs work...

(To be clear, by 'it' here I meant np.random.set_seed(), not the whole np.random API. Probably. And by 'deprecate' I mean 'whine loudly in some fashion when people use it', not 'rip out in a few releases'. I think.) -n

Nathaniel Smith wrote:

On Tue, Mar 12, 2013 at 9:25 PM, Nathaniel Smith wrote:

...

On Mon, Mar 11, 2013 at 9:46 AM, Robert Kern wrote:

...

On Sun, Mar 10, 2013 at 6:12 PM, Siu Kwan Lam wrote:

...

My suggestion to overcome (1) and (2) is to allow the user to select between the two implementations (and possibly different algorithms in the future). If user does not provide a choice, we use the MT19937-32 by default.

numpy.random.set_state("MT19937_64", …) # choose the 64-bit implementation

Most likely, the different PRNGs should be different subclasses of RandomState. The module-level convenience API should probably be left alone. If you need to control the PRNG that you are using, you really need to be passing around a RandomState instance and not relying on reseeding the shared global instance.

+1

...

Aside: I really wish we hadn't exposed `set_state()` in the module API. It's an attractive nuisance.

And our own test suite is a serious offender in this regard, we have tests that fail if you run the test suite in a non-default order... https://github.com/numpy/numpy/issues/347

I wonder if we dare deprecate it? The whole idea of a global random state is just a bad one, like every other sort of global shared state. But it's one that's deeply baked into a lot of scientific programmers expectations about how APIs work...

(To be clear, by 'it' here I meant np.random.set_seed(), not the whole np.random API. Probably. And by 'deprecate' I mean 'whine loudly in some fashion when people use it', not 'rip out in a few releases'. I think.)

-n

What do you mean that the idea of global shared state is a bad one? How would you prefer the API to look? An alternative is a stateless rng, where you have to pass it it's state on each invocation, which it would update and return. I hope you're not advocating that.

I guess I talked to you about 100 years ago about sharing state between numpy rng and code I have in c++ that wraps boost::random. So is there a C-api for this RandomState object I could use to call from c++? Maybe I could do something with that. The c++ code could invoke via the python api, but that might be slower. I'm just rambling here, I'd have to see the API to get some ideas.

Neal Becker wrote:

Neal Becker wrote:

...

I guess I talked to you about 100 years ago about sharing state between numpy rng and code I have in c++ that wraps boost::random. So is there a C-api for this RandomState object I could use to call from c++? Maybe I could do something with that.

The c++ code could invoke via the python api, but that might be slower. I'm just rambling here, I'd have to see the API to get some ideas.

I think if I could just grab a long int from the underlying mersenne twister, through some c api?

Well, this at least appears to work - probably not the most efficient approach - calls the RandomState object via the python interface to get 4 bytes at a time: int test1 (bp::object & rs) { bp::str bytes = call_methodbp::str (rs.ptr(), "bytes", 4); // get 4 bytes return *reinterpret_cast (PyString_AS_STRING (bytes.ptr())); } BOOST_PYTHON_MODULE (numpy_rand) { boost::numpy::initialize(); def ("test1", &test1); }

Robert Kern wrote:

On Wed, Mar 13, 2013 at 12:16 AM, Neal Becker wrote:

...

I guess I talked to you about 100 years ago about sharing state between numpy rng and code I have in c++ that wraps boost::random. So is there a C-api for this RandomState object I could use to call from c++? Maybe I could do something with that.

There is not one currently. Cython has provisions for sharing such low-level access to other Cython extensions, but I'm not sure how well it works for exporting data pointers and function pointers to general C/++ code. We could probably package the necessities into a struct and export a pointer to it via a PyCapsule.

I did find a way to do this, and the results are good enough. Timing is quite comparable to my pure c++ implementation. I used rk_ulong from mtrand.so. I also tried using rk_fill, but it was a bit slower. The boost::python c++ code is attached, for posterity.

On Thu, Mar 14, 2013 at 11:00 AM, Neal Becker wrote:

Robert Kern wrote:

...

On Wed, Mar 13, 2013 at 12:16 AM, Neal Becker wrote:

...

I guess I talked to you about 100 years ago about sharing state between numpy rng and code I have in c++ that wraps boost::random. So is there a C-api for this RandomState object I could use to call from c++? Maybe I could do something with that.

There is not one currently. Cython has provisions for sharing such low-level access to other Cython extensions, but I'm not sure how well it works for exporting data pointers and function pointers to general C/++ code. We could probably package the necessities into a struct and export a pointer to it via a PyCapsule.

One thing this code doesn't do: it requires construction of the wrapper class passing in a RandomState object. It doesn't verify you actually gave it a RandomState object. It's hard to do that. The problem as I see it is to perform this check, I need the RandomStateType object, which unfortunately mtrand.so does not export.

The only way to do it is in c++ code:

1. import numpy.random 2. get RandomState class 3. call it to create RandomState instance 4. get the ob_type pointer.

Pretty ugly:

object mod = object (handle<> (borrowed((PyImport_ImportModule("numpy.random"))))); object rs_obj = mod.attr("RandomState"); object rs_inst = call<object> (rs_obj.ptr(), 0); RandomStateTypeObj = rs_inst.ptr()->ob_type;

PyObject_IsInstance() should be sufficient. http://docs.python.org/2/c-api/object.html#PyObject_IsInstance -- Robert Kern

On Tue, Mar 12, 2013 at 7:10 PM, Robert Kern wrote:

On Tue, Mar 12, 2013 at 10:38 PM, Neal Becker wrote:

...

Nathaniel Smith wrote:

...

On Tue, Mar 12, 2013 at 9:25 PM, Nathaniel Smith wrote:

...

On Mon, Mar 11, 2013 at 9:46 AM, Robert Kern wrote:

...

On Sun, Mar 10, 2013 at 6:12 PM, Siu Kwan Lam wrote:

...

My suggestion to overcome (1) and (2) is to allow the user to select between the two implementations (and possibly different algorithms in the future). If user does not provide a choice, we use the MT19937-32 by default.

numpy.random.set_state("MT19937_64", …) # choose the 64-bit implementation

Most likely, the different PRNGs should be different subclasses of RandomState. The module-level convenience API should probably be left alone. If you need to control the PRNG that you are using, you really need to be passing around a RandomState instance and not relying on reseeding the shared global instance.

+1

...

Aside: I really wish we hadn't exposed `set_state()` in the module API. It's an attractive nuisance.

And our own test suite is a serious offender in this regard, we have tests that fail if you run the test suite in a non-default order... https://github.com/numpy/numpy/issues/347

I wonder if we dare deprecate it? The whole idea of a global random state is just a bad one, like every other sort of global shared state. But it's one that's deeply baked into a lot of scientific programmers expectations about how APIs work...

(To be clear, by 'it' here I meant np.random.set_seed(), not the whole np.random API. Probably. And by 'deprecate' I mean 'whine loudly in some fashion when people use it', not 'rip out in a few releases'. I think.)

-n

What do you mean that the idea of global shared state is a bad one?

The words "global shared state" drives fear into the hearts of experienced programmers everywhere, whatever the context. :-) It's rarely a *good* idea.

...

How would you prefer the API to look?

There are two current APIs:

1. Instantiate RandomState and call it's methods 2. Just call the functions in numpy.random

The latter has a shared global state. In fact, all of those "functions" are just references to the methods on a shared global RandomState instance.

We advocate using the former API. Note that it already exists. It was the recommended API from day one. No one is recommending adding a new API.

I never saw much advertising for the RandomState api, and until recently wasn't sure why using the global random state function np.random.norm, ... should be a bad idea. Learning by example, and seeing almost all examples using the global state, is not exactly conducive to figuring out that there is an issue. All of scipy.stats.distribution random numbers are using the global random state. (I guess I should open a ticket.) Josef

...

An alternative is a stateless rng, where you have to pass it it's state on each invocation, which it would update and return. I hope you're not advocating that.

No. This is a place where OOP solved the problem neatly.

-- Robert Kern _______________________________________________ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion