Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Forms of propensity calculations #35

Open
wants to merge 7 commits into
base: master
Choose a base branch
from
Open

Forms of propensity calculations #35

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
7 commits
Select commit Hold shift + click to select a range
8985e96
Add optional 'forms' input argument.
heejochoi Feb 7, 2019
df877d3
Add test for optional 'forms' input argument.
heejochoi Feb 7, 2019
5a60cd9
Switch propensity computation based on reaction form.
heejochoi Feb 8, 2019
62d037a
Update docstrings.
heejochoi Feb 8, 2019
de698ef
Restructure 'rates' into a list of variable length lists.
heejochoi Feb 9, 2019
cf462cc
Change constant in fraction saturation computation to double.
heejochoi Feb 9, 2019
cc2c522
Update test_forms.
heejochoi Feb 9, 2019
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
25 changes: 18 additions & 7 deletions arrow/arrow.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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):
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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?

'''
This invokes the Obsidian C module (see obsidianmodule.c) with the
stoichiometry, reaction rates and a variety of derived values. Once constructed,
this can be invoked by calling `evolve` with a duration and initial state, since
the stoichiometry will be shared among all calls.
stoichiometry, reaction rates, the form of each reaction (optional)
and a variety of derived values. Once constructed, this can be invoked
by calling `evolve` with a duration and initial state, since the
stoichiometry will be shared among all calls.

There are four derived values, each of which is a list of variable length
lists. In order to pass this into C, these nested lists are flattened and two
Expand All @@ -101,7 +102,14 @@ def __init__(self, stoichiometry, rates, random_seed=0):
'''

self.stoichiometry = stoichiometry
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)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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

  1. I think there is a more NumPy-flavored solution to converting a 1D array to a 2D array - see np.atleast_2d.
  2. 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?


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

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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.


self.random_seed = random_seed

Expand All @@ -117,7 +125,10 @@ def __init__(self, stoichiometry, rates, random_seed=0):
self.obsidian = obsidian.obsidian(
self.random_seed,
self.stoichiometry,
self.rates,
self.rates_flat,
self.rates_lengths,
self.rates_indexes,
self.forms,
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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.

self.reactants_lengths,
self.reactants_indexes,
self.reactants_flat,
Expand All @@ -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])
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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

for event in events:
occurrences[event] += 1

Expand Down
73 changes: 57 additions & 16 deletions arrow/obsidian.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -24,6 +24,12 @@ choose(long n, long k) {
return combinations;
}

// Compute the fractional saturation
double
fraction_saturation(long s, double k) {
return ((double) s / (s + k));
}

// Perform the Gillespie algorithm with the given stoichiometry, reaction rates and
// initial state for the provided duration.

Expand All @@ -34,11 +40,14 @@ choose(long n, long k) {
// contains all the information about the reactions this system performs. Each
// row is a reaction, and each column is a substrate, so every entry is the
// change in counts of each substrate when the given reaction is performed.
// * rates: an array of length `reactions_count` that encodes the base rate for
// each reaction. The actual propensity is further dependent on the counts for
// each reactant in the reaction.
// * rates_flat: an array of length `reactions_count` (or greater if the reaction
// is not of the default form 0) that encodes the base rate for each reaction.
// The actual propensity is further dependent on the counts for each reactant
// in the reaction.
// * forms: an array of length `reactions_count` that indicates the type of
// propensity computation that is associated with each reaction.
// * duration: How long to run the simulation for.
// * state: An array of length `substrates_count` that contains the inital count
// * state: An array of length `substrates_count` that contains the initial count
// for each substrate. The outcome of the algorithm will involve an array of
// the same size that signifies the counts of each substrate after all the
// reactions are performed.
Expand Down Expand Up @@ -83,7 +92,10 @@ evolve(MTState *random_state,
int reactions_count,
int substrates_count,
long *stoichiometry,
double *rates,
double *rates_flat,
long *rates_lengths,
long *rates_indexes,
long *forms,

long *reactants_lengths,
long *reactants_indexes,
Expand Down Expand Up @@ -158,18 +170,47 @@ evolve(MTState *random_state,
// First update all propensities that were affected by the previous event
for (up = 0; up < update_length; up++) {

// Find which reaction to update and initialize the propensity for that reaction
// with the rate for that reaction
// Find which reaction to update
reaction = update[up];
propensities[reaction] = rates[reaction];

// Go through each reactant and calculate its contribution to the propensity
// based on the counts of the corresponding substrate, which is then multiplied
// by the reaction's original rate and the contributions from other reactants
for (reactant = 0; reactant < reactants_lengths[reaction]; reactant++) {
index = reactants_indexes[reaction] + reactant;
count = outcome[reactants[index]];
propensities[reaction] *= choose(count, reactions[index]);

// Compute propensity depending on the reaction's form
switch(forms[reaction]) {

case 0: // Standard Gillespie
// Initialize the propensity for reaction with its rate constant
propensities[reaction] = rates_flat[reaction];

// Go through each reactant and calculate its contribution to the
// propensity based on the counts of the corresponding substrate,
// which is then multiplied by the reaction's original rate and the
// contributions from other reactants
for (reactant = 0; reactant < reactants_lengths[reaction]; reactant++) {
index = reactants_indexes[reaction] + reactant;
count = outcome[reactants[index]];
propensities[reaction] *= choose(count, reactions[index]);
}

break;

case 1: // Michaelis-Menten
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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.

// Initialize the propensity for reaction with its maximal propensity
index = rates_indexes[reaction];
propensities[reaction] = rates_flat[index];

// Go through each reactant and calculate its contribution to the
// propensity based on the counts of the corresponding substrates
// relative to their Michaelis-Menten constant.
for (reactant = 0; reactant < reactants_lengths[reaction]; reactant++) {
index = reactants_indexes[reaction] + reactant;
count = outcome[reactants[index]];
propensities[reaction] *= fraction_saturation(count,
rates_flat[rates_indexes[reaction] + reactant + 1]);
}

break;

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

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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

}
}

Expand Down
5 changes: 4 additions & 1 deletion arrow/obsidian.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -19,7 +19,10 @@ evolve(MTState *random_state,
int reactions_count,
int substrates_count,
long *stoichiometry,
double *rates,
double *rates_flat,
long *rates_lengths,
long *rates_indexes,
long *forms,

long *reactants_lengths,
long *reactants_indexes,
Expand Down
55 changes: 44 additions & 11 deletions arrow/obsidianmodule.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -39,7 +39,10 @@ typedef struct {
int reactions_count;
int substrates_count;
long *stoichiometry;
double *rates;
double *rates_flat;
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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.

long *rates_lengths;
long *rates_indexes;
long *forms;

long *reactants_lengths;
long *reactants_indexes;
Expand All @@ -66,7 +69,10 @@ newObsidianObject(int random_seed,
int reactions_count,
int substrates_count,
long *stoichiometry,
double *rates,
double *rates_flat,
long *rates_lengths,
long *rates_indexes,
long *forms,

long *reactants_lengths,
long *reactants_indexes,
Expand Down Expand Up @@ -103,7 +109,10 @@ newObsidianObject(int random_seed,
self->reactions_count = reactions_count;
self->substrates_count = substrates_count;
self->stoichiometry = stoichiometry;
self->rates = rates;
self->rates_flat = rates_flat;
self->rates_lengths = rates_lengths;
self->rates_indexes = rates_indexes;
self->forms = forms;

self->reactants_lengths = reactants_lengths;
self->reactants_indexes = reactants_indexes;
Expand Down Expand Up @@ -138,7 +147,7 @@ Obsidian_demo(ObsidianObject *self, PyObject *args)
if (!PyArg_ParseTuple(args, ":demo"))
return NULL;

print_array(self->rates, self->reactions_count);
print_array(self->rates_flat, self->reactions_count);

Py_INCREF(Py_None);
return Py_None;
Expand Down Expand Up @@ -187,7 +196,10 @@ Obsidian_evolve(ObsidianObject *self, PyObject *args)
self->reactions_count,
self->substrates_count,
self->stoichiometry,
self->rates,
self->rates_flat,
self->rates_lengths,
self->rates_indexes,
self->forms,
self->reactants_lengths,
self->reactants_indexes,
self->reactants,
Expand Down Expand Up @@ -364,7 +376,10 @@ _invoke_obsidian(PyObject *self, PyObject *args) {

int random_seed;
PyObject *stoichiometry_obj,
*rates_obj,
*rates_flat_obj,
*rates_lengths_obj,
*rates_indexes_obj,
*forms_obj,

*reactants_lengths_obj,
*reactants_indexes_obj,
Expand All @@ -380,10 +395,13 @@ _invoke_obsidian(PyObject *self, PyObject *args) {
*substrates_obj;

if (!PyArg_ParseTuple(args,
"iOOOOOOOOOOOO",
"iOOOOOOOOOOOOOOO",
&random_seed,
&stoichiometry_obj,
&rates_obj,
&rates_flat_obj,
&rates_lengths_obj,
&rates_indexes_obj,
&forms_obj,

&reactants_lengths_obj,
&reactants_indexes_obj,
Expand All @@ -406,8 +424,16 @@ _invoke_obsidian(PyObject *self, PyObject *args) {
long *stoichiometry = (long *) PyArray_DATA(stoichiometry_array);

// Import the rates for each reaction
PyObject *rates_array = array_for(rates_obj, NPY_DOUBLE);
double *rates = (double *) PyArray_DATA(rates_array);
PyObject *rates_flat_array = array_for(rates_flat_obj, NPY_DOUBLE);
double *rates_flat = (double *) PyArray_DATA(rates_flat_array);
PyObject *rates_lengths_array = array_for(rates_lengths_obj, NPY_INT64);
long *rates_lengths = (long *) PyArray_DATA(rates_lengths_array);
PyObject *rates_indexes_array = array_for(rates_indexes_obj, NPY_INT64);
long *rates_indexes = (long *) PyArray_DATA(rates_indexes_array);

// Import the forms for each reaction
PyObject *forms_array = array_for(forms_obj, NPY_INT64);
long *forms = (long *) PyArray_DATA(forms_array);

// Pull all the precalculated nested arrays
PyObject *reactants_lengths_array = array_for(reactants_lengths_obj, NPY_INT64);
Expand Down Expand Up @@ -438,7 +464,10 @@ _invoke_obsidian(PyObject *self, PyObject *args) {
reactions_count,
substrates_count,
stoichiometry,
rates,
rates_flat,
rates_lengths,
rates_indexes,
forms,

reactants_lengths,
reactants_indexes,
Expand All @@ -459,6 +488,10 @@ _invoke_obsidian(PyObject *self, PyObject *args) {
}

// Clean up all the PyObject * references
Py_XDECREF(rates_flat_array);
Py_XDECREF(rates_lengths_array);
Py_XDECREF(rates_indexes_array);

Py_XDECREF(reactants_lengths_array);
Py_XDECREF(reactants_indexes_array);
Py_XDECREF(reactants_array);
Expand Down
19 changes: 19 additions & 0 deletions arrow/test/test_arrow.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -173,6 +173,25 @@ def test_compare_runtime():
print('reference time elapsed: {}'.format(reference_end - reference_start))
print('obsidian time elapsed: {}'.format(obsidian_end - obsidian_start))

def test_forms():
stoichiometric_matrix = np.array([
[-1, -1, 1, 0, 0, 0],
[0, 0, 0, -1, -1, 1]], np.int64)

rates = np.array([
[10, 1, 1],
[1, 1, 1]], np.float64)
forms = np.array([1, 1])

arrow = StochasticSystem(stoichiometric_matrix, rates, forms)
result = arrow.evolve(1., np.array([50, 50, 0, 20, 20, 0], np.int64))

print('steps: {}'.format(result['steps']))
print('events: {}'.format(result['events']))
print('occurrences: {}'.format(result['occurrences']))
print('outcome: {}'.format(result['outcome']))


if __name__ == '__main__':
parser = argparse.ArgumentParser()
parser.add_argument('--plot', action='store_true')
Expand Down