Merge pull request #49 from CovertLab/fix-48

Avoid selection of 0 propensity reactions and negative counts

Author: tahorst

README.md

@@ -104,8 +104,25 @@ There are more command line features in test_arrow:
 
     > python -m arrow.test.test_arrow --time
 
+More examples:
+
+    > python -m arrow.test.test_hang
+
+    > pytest -m arrow/test/test_arrow.py
+
+    > pytest -k flagella
+
 ## Changelog
 
+### Version 0.5.0
+
+* Add the arrow_hang unit test which catches a nasty edge-case (Issue #48),
+  fix the bug, and make the code more robust to some other potential bugs.
+
+### Version 0.4.4
+
+* Can pickle StochasticSystem instances.
+
 ### Version 0.3.0
 
 * Introduced backwards-incompatible API change for supplying rates at `evolve()` time rather than `__init__()` for `StochasticSystem`.

arrow/obsidian.c

@@ -175,9 +175,16 @@ evolve_result evolve(Info *info, double duration, int64_t *state, double *rates)
     // Find the total for all propensities
     total = 0.0;
     for (reaction = 0; reaction < reactions_count; reaction++) {
+      if (propensities[reaction] < 0) {
+        status = 3; // a negative propensity
+      }
       total += propensities[reaction];
     }
 
+    if (status > 0) {
+      break;
+    }
+
     if (isnan(total)) {
       printf("failed simulation: total propensity is NaN\n");
       int max_reaction = 0;
@@ -191,6 +198,7 @@ evolve_result evolve(Info *info, double duration, int64_t *state, double *rates)
       interval = 0.0;
       choice = -1;
       status = 1; // overflow
+      break;
     }
 
     // If the total is zero, then we have no more reactions to perform and can exit
@@ -207,22 +215,29 @@ evolve_result evolve(Info *info, double duration, int64_t *state, double *rates)
       // `interval` from an exponential distribution.
       interval = sample_exponential(random_state, total);
       point = sample_uniform(random_state) * total;
+      if (point > total) {
+        // If roundoff made point > total, the `progress` loop would go past the
+        // end of the array.
+        point = total;
+      }
+
+      // If we have surpassed the provided duration we can exit now
+      if (now + interval > duration) {
+        break;
+      }
 
       // Based on the random sample, find the event that it corresponds to by
       // iterating through the propensities until we surpass our sampled value
       choice = 0;
       progress = 0.0;
 
-      while (progress + propensities[choice] < point) {
+      // Note: Even if `point` happens to be 0, this needs to skip 0 propensity
+      // choices to avoid computing negative counts.
+      while (progress + propensities[choice] < point || propensities[choice] == 0) {
         progress += propensities[choice];
         choice += 1;
       }
 
-      // If we have surpassed the provided duration we can exit now
-      if (choice == -1 || (now + interval) > duration) {
-        break;
-      }
-
       // Increase time by the interval sampled above
       now += interval;
 
@@ -237,6 +252,14 @@ evolve_result evolve(Info *info, double duration, int64_t *state, double *rates)
         index = substrates_indexes[choice] + involve;
         adjustment = stoichiometry[choice * substrates_count + substrates[index]];
         outcome[substrates[index]] += adjustment;
+
+        if (outcome[substrates[index]] < 0) {
+          status = 2; // negative counts
+        }
+      }
+
+      if (status > 0) {
+        break;
       }
 
       // Find which propensities depend on this reaction and therefore need to be

arrow/test/complex-counts.npy

@@ -0,0 +1,56 @@
+# A test case for Issue #48.
+#
+# This code should reproduce the error: the program hangs after it prints "7100"
+# and eventually runs out of memory. You'll have to Ctrl+Z to stop it.
+# Ctrl+C won't work.
+#
+# It hangs only with certain seeds and numbers of molecules. The system can
+# evolve with the same number of molecule counts for 7179 iterations before it
+# hangs. Adding 1 to all of the molecules causes it to hang at an earlier
+# iteration.
+#
+# TODO: Debug this. Is it caused by a Gillespie algorithm blowup (see below),
+# integer overflow Undefined Behavior in C, or something else?
+#
+# The Gillespie algorithm is prone to explode [symptom?] under certain
+# conditions if the exponent term in the choice calculation is too large.
+#
+# The workaround is to find the problematic reaction and decompose the
+# stoichiometry into an equivalent problem with more steps.
+#
+# The flagella example had something like 170 identical subunits which caused
+# the problem. Breaking it into 2+ equivalent reactions fixed it.
+#
+# It'd be good for the Arrow code to catch this problem when/before it happens
+# and at least identify which reactions are problematic.
+
+import os
+
+from arrow import StochasticSystem
+import numpy as np
+
+
+def np_load(filename):
+    filepath = os.path.join(os.path.dirname(__file__), filename)
+    return np.load(filepath)
+
+
+def test_hang():
+    # TODO: Use a pytest plug-in to timeout after some threshold.
+
+    seed = 807952948
+
+    stoich = np_load('stoich.npy')
+    mol = np_load('complex-counts.npy')
+    rates = np_load('rates.npy')
+
+    system = StochasticSystem(stoich, random_seed=seed)
+    for i in range(7300):
+        if i % 100 == 0:
+            print(i)
+
+        result = system.evolve(1, mol, rates)
+
+
+if __name__ == '__main__':
+    test_hang()

arrow/test/rates.npy

@@ -1,12 +1,11 @@
 import os
-# from glob import glob
 import setuptools  # used indirectly for bdist_wheel cmd and long_description_content_type
 from distutils.core import setup
 from distutils.extension import Extension
 import numpy.distutils.misc_util
 
 with open("README.md", 'r') as readme:
-	long_description = readme.read()
+    long_description = readme.read()
 
 current_dir = os.getcwd()
 arrow_dir = os.path.join(current_dir, 'arrow')
@@ -23,38 +22,38 @@
 ext = '.pyx' if USE_CYTHON else '.c'
 
 cython_extensions = [
-	Extension('arrow.arrowhead',
-			  sources=['arrow/mersenne.c', 'arrow/obsidian.c', 'arrow/arrowhead'+ext,],
-			  include_dirs=['arrow'] + numpy.distutils.misc_util.get_numpy_include_dirs(),
-			  define_macros=[('NPY_NO_DEPRECATED_API', 'NPY_1_7_API_VERSION')],
-			  )]
+    Extension('arrow.arrowhead',
+              sources=['arrow/mersenne.c', 'arrow/obsidian.c', 'arrow/arrowhead'+ext,],
+              include_dirs=['arrow'] + numpy.distutils.misc_util.get_numpy_include_dirs(),
+              define_macros=[('NPY_NO_DEPRECATED_API', 'NPY_1_7_API_VERSION')],
+              )]
 
 if USE_CYTHON:
-	from Cython.Build import cythonize
-	cython_extensions = cythonize(
-		cython_extensions,
-		include_path=['arrow'],
-		annotate=True,  # to get an HTML code listing
-	  	)
+    from Cython.Build import cythonize
+    cython_extensions = cythonize(
+        cython_extensions,
+        include_path=['arrow'],
+        annotate=True,  # to get an HTML code listing
+    )
 
 setup(
-	name='stochastic-arrow',
-	version='0.4.4',
-	packages=['arrow'],
-	author='Ryan Spangler, John Mason, Jerry Morrison',
-	author_email='spanglry@stanford.edu',
-	url='https://github.com/CovertLab/arrow',
-	license='MIT',
-	include_dirs=include,
-	ext_modules=cython_extensions,
-	long_description=long_description,
-	long_description_content_type='text/markdown',
-	requires=['numpy (>=1.14)', 'six'],
-	classifiers=[
-		'Development Status :: 3 - Alpha',
-		'License :: OSI Approved :: MIT License',
-		'Programming Language :: Python',
-		'Programming Language :: Python :: 2.7',
-		'Programming Language :: Python :: 3',
-		'Topic :: Scientific/Engineering',
-	])
+    name='stochastic-arrow',
+    version='0.5.0',
+    packages=['arrow'],
+    author='Ryan Spangler, John Mason, Jerry Morrison, Chris Skalnik, Travis Ahn-Horst',
+    author_email='ryan.spangler@gmail.com',
+    url='https://github.com/CovertLab/arrow',
+    license='MIT',
+    include_dirs=include,
+    ext_modules=cython_extensions,
+    long_description=long_description,
+    long_description_content_type='text/markdown',
+    requires=['numpy (>=1.14)', 'six'],
+    classifiers=[
+        'Development Status :: 3 - Alpha',
+        'License :: OSI Approved :: MIT License',
+        'Programming Language :: Python',
+        'Programming Language :: Python :: 2.7',
+        'Programming Language :: Python :: 3',
+        'Topic :: Scientific/Engineering',
+    ])