Forms of propensity calculations #35

heejochoi · 2019-02-13T00:56:04Z

This PR aims to solve issue #34 by facilitating alternative propensity calculations:

Adds optional input argument called forms, which mediates using different functions to calculate the reaction's propensity. If forms is not provided as an input argument, the propensities of all reactions in the system are computed according to standard Gillespie -- ie. defaults to current behavior.
Allows the existing rates input argument to be optionally a list (or array) of nested lists (or arrays) of variable lengths to accommodate for propensity calculations that involve more than one constant.
Adds test for forms.

prismofeverything

Hey @heejochoi, looks good! I couldn't find any issues with this, I ran the code and the tests all pass. I have a couple comments but mostly just planning for the future.

prismofeverything · 2019-02-13T21:42:03Z

arrow/arrow.py

-        self.rates = rates
+        if rates.ndim == 1:
+            rates = [[rate] for rate in rates.astype(np.float64)]
+        self.rates_flat, self.rates_lengths, self.rates_indexes = flat_indexes(rates)


Good! This means our rates are always in the same form, even if there is just a single rate per reaction.

Apart from not being clear why 'rates' would be something other than 1D,

I think there is a more NumPy-flavored solution to converting a 1D array to a 2D array - see np.atleast_2d.

Why are you casting to floats in here? If it's not 1D, is it OK if they are not floats? Do they have to be floats at all, and if so, could that be handled at a different layer?

prismofeverything · 2019-02-13T21:44:01Z

arrow/arrow.py

+        if forms is not None:
+            self.forms = forms
+        else:
+            self.forms = np.zeros(stoichiometry.shape[0], dtype=np.int64)


Again, good. Leave all the optional functionality in the python and pass one consistent form to the C. Thanks for this.

I think it's a bit more elegant to write code like

if forms is None: forms = # some sort of default self.forms = forms

which separates out the input parsing from the attribute assignment. However this is more my preference than a standard, if I'm not mistaken. This is a common enough case that you'd think there would be a standard but I've never seen mention of one.

prismofeverything · 2019-02-13T21:45:15Z

arrow/arrow.py

+            self.rates_flat,
+            self.rates_lengths,
+            self.rates_indexes,
+            self.forms,


More args! Let's talk about what it looks like to make structs for the nested lists.... may be more trouble than it is worth, but could also clarify some things. It is a tradeoff we should evaluate.

prismofeverything · 2019-02-13T21:47:31Z

arrow/obsidian.c

+
+          break;
+
+        case 1: // Michaelis-Menten


Right, the form switch. We may want to talk about what it looks like to

Add more forms.

Use user provided forms.

Move the essence of each form into its own function?

There will get to be a point where this switch gets unwieldy, but it serves for now.

prismofeverything · 2019-02-13T21:48:30Z

arrow/obsidianmodule.c

@@ -39,7 +39,10 @@ typedef struct {
  int reactions_count;
  int substrates_count;
  long *stoichiometry;
-  double *rates;
+  double *rates_flat;


Love how many places you have to add this.... Looking forward to evaluating @1fish2's proposal of using cython for the glue (module) code here.

jmason42

Coming back to this after a while is interesting, as it highlights the places where the code isn't very clear. I've added some comments to that effect.

From what I can glean, you've implemented a pseudo-convenience kinetics form as an alternative to the standard stochastic form. I don't want to throw up any roadblocks here but I do want to make it clear that anyone who uses this form is playing with fire, for several reasons:

Michaelis-Menten kinetics, and its derivatives, are all based on concentrations, while stochastic simulation algorithms utilize molecule counts. You'll need to be very careful about converting between kinetic rate constants and their stochastic equivalents, including your new saturation constants (i.e. Michaelis-Menten constants).
This pseudo-convenience kinetics approach will break in a case where a reactant has a stoichiometry of 2 if the number of molecules goes to 1.
Using a saturation function in a stochastic simulation is on shaky grounds compared to the original form. I don't think it's totally out of the question - you're bridging between two modeling approaches and that often requires trying things before worrying about correctness.
You're not handling any type of competition, including the competitive effects between the forward and reverse reactants for an enzyme catalyzing a reversible reaction. So, be very apprehensive about introducing any reversible reactions. Otherwise you'll be double-counting enzyme capacity.

jmason42 · 2019-02-14T19:35:51Z

arrow/arrow.py

@@ -80,12 +80,13 @@ class StochasticSystem(object):
    (and zero everywhere else). 
    '''

-    def __init__(self, stoichiometry, rates, random_seed=0):
+    def __init__(self, stoichiometry, rates, forms=None, random_seed=0):


This needs a doc string, as I've forgotten what 'rates' is and don't know where to look. :)

Are these stochastic rate constants, perhaps? How do they play with the new 'forms'? Are 'forms' like kinetic rate laws sans the outermost proportionality constant?

jmason42 · 2019-02-14T19:41:28Z

arrow/arrow.py

-        self.rates = rates
+        if rates.ndim == 1:
+            rates = [[rate] for rate in rates.astype(np.float64)]
+        self.rates_flat, self.rates_lengths, self.rates_indexes = flat_indexes(rates)


Apart from not being clear why 'rates' would be something other than 1D,

I think there is a more NumPy-flavored solution to converting a 1D array to a 2D array - see np.atleast_2d.

Why are you casting to floats in here? If it's not 1D, is it OK if they are not floats? Do they have to be floats at all, and if so, could that be handled at a different layer?

jmason42 · 2019-02-14T19:44:30Z

arrow/arrow.py

+        if forms is not None:
+            self.forms = forms
+        else:
+            self.forms = np.zeros(stoichiometry.shape[0], dtype=np.int64)


I think it's a bit more elegant to write code like

if forms is None: forms = # some sort of default self.forms = forms

which separates out the input parsing from the attribute assignment. However this is more my preference than a standard, if I'm not mistaken. This is a common enough case that you'd think there would be a standard but I've never seen mention of one.

jmason42 · 2019-02-14T19:45:45Z

arrow/arrow.py

@@ -131,7 +142,7 @@ def __init__(self, stoichiometry, rates, random_seed=0):

    def evolve(self, duration, state):
        steps, time, events, outcome = self.obsidian.evolve(duration, state)
-        occurrences = np.zeros(len(self.rates))
+        occurrences = np.zeros(self.stoichiometry.shape[0])


I think self.stoichiometry.shape[0] is used in a couple places and should probably be factored out for convenience/consistency.

jmason42 · 2019-02-14T19:48:59Z

arrow/obsidian.c

+          break;
+
+        default:
+          printf("arrow.obsidian.evolve - unexpected form: %ld", forms[reaction]);


This probably needs to throw an error rather than just print a message.

heejochoi added 7 commits February 7, 2019 11:19

Add optional 'forms' input argument.

8985e96

Add test for optional 'forms' input argument.

df877d3

Switch propensity computation based on reaction form.

5a60cd9

Update docstrings.

62d037a

Restructure 'rates' into a list of variable length lists.

de698ef

Change constant in fraction saturation computation to double.

cf462cc

Update test_forms.

cc2c522

prismofeverything approved these changes Feb 13, 2019

View reviewed changes

jmason42 reviewed Feb 14, 2019

View reviewed changes

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Forms of propensity calculations #35

Forms of propensity calculations #35

heejochoi commented Feb 13, 2019

prismofeverything left a comment

prismofeverything Feb 13, 2019

jmason42 Feb 14, 2019

prismofeverything Feb 13, 2019

jmason42 Feb 14, 2019

prismofeverything Feb 13, 2019

prismofeverything Feb 13, 2019

prismofeverything Feb 13, 2019

jmason42 left a comment

jmason42 Feb 14, 2019

jmason42 Feb 14, 2019

jmason42 Feb 14, 2019

jmason42 Feb 14, 2019

jmason42 Feb 14, 2019

jmason42 Feb 14, 2019

Forms of propensity calculations #35

Are you sure you want to change the base?

Forms of propensity calculations #35

Conversation

heejochoi commented Feb 13, 2019

prismofeverything left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

jmason42 left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment