Implement Gillespie in C and call from python

[prismofeverything]

This PR dispenses struggling with trying to make python fast and appeals directly to C, the source of all things.

The python interface is the same (though it returns a dictionary instead of a tuple as the values returned were growing), but instead of using the now reference python implementation the arguments are translated from python into C arrays and passed to a pure C Gillespie implementation. Comparing the two:

reference time elapsed: 2.04260396957
obsidian time elapsed: 0.0289080142975

This is in the range of two orders of magnitude, or over 70x faster. Luckily for us, the Gillespie algorithm can never run too fast, so some day I hope to see this realized: https://www.ncbi.nlm.nih.gov/pubmed/19162916 (anyone know who the FPGA people are on campus? I have already asked several of you). Until then, we have C. Which will suffice for now.

I go into great detail about the design and algorithm throughout the code so I will avoid that here. The main entry point and translation from python occurs in obsidianmodule.c and the Gillespie implementation is in obsidian.c with a header at obsidian.h, in the evolve function. Any and all questions/suggestions/thoughts/revelations welcome, thanks!

[prismofeverything]

Thanks @1fish2 and @jmason42! Savvy feedback as always, much appreciated.

Easy, high quality fix: Call r.rand(1000) to get (e.g.) 1000 random numbers where r = numpy.random.RandomState(seed) with the caller's seed.

Right, I think these can be generated and passed in, though it means adding another argument to the function. I wonder also if just the seed can be passed in and some deterministic RNG could be used directly. It needs two random numbers per step, so it seems we can't get around having to generate them on the fly if we exceed the number of pregenerated ones. I'll look into a deterministic RNG in C.... usually you get the random number and the next seed, which can be stored and passed around.

Let's require new code to be compatible with Python 3. Is this code already good there? See Porting Extension Modules to Python 3.

Ah yes, I did start this before we decided to go all out on python 3. That said it doesn't look like the changes are that significant, though this could be the opportunity to let cython do the work as you say. I'll look into that at some point.

[1fish2]

Yes, pass in a np.RandomState, or pass in a seed and call NumPy, or pass in a seed and then call its own random number generator, e.g. crib NumPy's Mersenne Twister. But finding, linking in, and testing a portable, high-quality generator sounds like a rabbit hole.

[prismofeverything]

Okay, found a good BSD-licensed mersenne twister in C++ I was able to adapt to our purposes: https://github.com/cslarsen/mersenne-twister

The random state was static before so I added it to the obsidian object during initialization and added it as an argument to all the mersenne functions. Now each new obsidian object will have its own random state that it can track and keep independent of other states throughout the simulation. I also added a few functions, sample_uniform and sample_expontential to convert the unsigned ints into the doubles we need: https://github.com/CovertLab/arrow/commit/90dc3f58dfed2ba96bf79e5c0cd20f023d5bd93d

The interface has acquired a new argument to pass in the seed.

I ran it a few times and the results are entirely identical. I haven't tested it with the sim yet but I assume this holds true.

[prismofeverything]

Okay, handling malloc errors now, what do you think @1fish2? https://github.com/CovertLab/arrow/commit/e0017e8ecef0d390ae07a7c8f39c6bf7e6a5d9a7

I'm not sure if I need to free the other allocated items if it fails since that means the program is likely ending soon anyway, but thought I would be careful there. The way to signify a python call has failed is to return NULL from it, so this will escalate to an error in the module code if malloc fails here, which I think is the correct behavior. Open to feedback here of course.

[prismofeverything]

Beyond that I think I am going to leave the module/cython/python 3 rehabilitation for another PR and circle around to that when we are ready. Otherwise I think I've dealt with most of the feedback here, thanks again for the thorough review : ) It is good to get some eyes on it that can be critical about the C practices, all of my C knowledge comes from 2 years in a C++ shop with a highly idiosyncratic codebase (the guy had basically implemented a whole operating system inside of an activeX plugin that only ran in IE6) so good to hear other ways of doing things : )

[prismofeverything]

FYI, free(NULL) is defined to be a no-op by the API spec, so these NULL checks are not needed. No harm but clutter.

Ah good to know, I'll adjust that.

The printfs are a good idea. I'd add a bit more context since the printout will appear out of the blue, e.g. "arrow.obsidian.evolve(): failed to allocate memory (%d)", errno. A developer seeing an allocation failure probably won't need to track down the source code, but they did they'll need more context than evolve.

Good call, I'll add that in. Yeah just so long as someone can figure out what actually went wrong.