From 8985e96b42d6f6863c51a96a94318781782c6e09 Mon Sep 17 00:00:00 2001 From: Heejo Choi Date: Thu, 7 Feb 2019 11:19:22 -0800 Subject: [PATCH 1/7] Add optional 'forms' input argument. --- arrow/arrow.py | 8 +++++++- arrow/obsidian.c | 1 + arrow/obsidian.h | 1 + arrow/obsidianmodule.c | 13 ++++++++++++- 4 files changed, 21 insertions(+), 2 deletions(-) diff --git a/arrow/arrow.py b/arrow/arrow.py index a7d3183..5b9981c 100644 --- a/arrow/arrow.py +++ b/arrow/arrow.py @@ -80,7 +80,7 @@ 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 invokes the Obsidian C module (see obsidianmodule.c) with the stoichiometry, reaction rates and a variety of derived values. Once constructed, @@ -103,6 +103,11 @@ def __init__(self, stoichiometry, rates, random_seed=0): self.stoichiometry = stoichiometry self.rates = rates + if forms is not None: + self.forms = forms + else: + self.forms = np.zeros(stoichiometry.shape[0], dtype=np.int64) + self.random_seed = random_seed reactants, reactions, substrates = derive_reactants(stoichiometry) @@ -118,6 +123,7 @@ def __init__(self, stoichiometry, rates, random_seed=0): self.random_seed, self.stoichiometry, self.rates, + self.forms, self.reactants_lengths, self.reactants_indexes, self.reactants_flat, diff --git a/arrow/obsidian.c b/arrow/obsidian.c index ab77a33..880e50d 100644 --- a/arrow/obsidian.c +++ b/arrow/obsidian.c @@ -84,6 +84,7 @@ evolve(MTState *random_state, int substrates_count, long *stoichiometry, double *rates, + long *forms, long *reactants_lengths, long *reactants_indexes, diff --git a/arrow/obsidian.h b/arrow/obsidian.h index ad70f6a..f2f0ea4 100644 --- a/arrow/obsidian.h +++ b/arrow/obsidian.h @@ -20,6 +20,7 @@ evolve(MTState *random_state, int substrates_count, long *stoichiometry, double *rates, + long *forms, long *reactants_lengths, long *reactants_indexes, diff --git a/arrow/obsidianmodule.c b/arrow/obsidianmodule.c index 4e4ee1d..c20752c 100644 --- a/arrow/obsidianmodule.c +++ b/arrow/obsidianmodule.c @@ -40,6 +40,7 @@ typedef struct { int substrates_count; long *stoichiometry; double *rates; + long *forms; long *reactants_lengths; long *reactants_indexes; @@ -67,6 +68,7 @@ newObsidianObject(int random_seed, int substrates_count, long *stoichiometry, double *rates, + long *forms, long *reactants_lengths, long *reactants_indexes, @@ -104,6 +106,7 @@ newObsidianObject(int random_seed, self->substrates_count = substrates_count; self->stoichiometry = stoichiometry; self->rates = rates; + self->forms = forms; self->reactants_lengths = reactants_lengths; self->reactants_indexes = reactants_indexes; @@ -188,6 +191,7 @@ Obsidian_evolve(ObsidianObject *self, PyObject *args) self->substrates_count, self->stoichiometry, self->rates, + self->forms, self->reactants_lengths, self->reactants_indexes, self->reactants, @@ -365,6 +369,7 @@ _invoke_obsidian(PyObject *self, PyObject *args) { int random_seed; PyObject *stoichiometry_obj, *rates_obj, + *forms_obj, *reactants_lengths_obj, *reactants_indexes_obj, @@ -380,10 +385,11 @@ _invoke_obsidian(PyObject *self, PyObject *args) { *substrates_obj; if (!PyArg_ParseTuple(args, - "iOOOOOOOOOOOO", + "iOOOOOOOOOOOOO", &random_seed, &stoichiometry_obj, &rates_obj, + &forms_obj, &reactants_lengths_obj, &reactants_indexes_obj, @@ -409,6 +415,10 @@ _invoke_obsidian(PyObject *self, PyObject *args) { PyObject *rates_array = array_for(rates_obj, NPY_DOUBLE); double *rates = (double *) PyArray_DATA(rates_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); long *reactants_lengths = (long *) PyArray_DATA(reactants_lengths_array); @@ -439,6 +449,7 @@ _invoke_obsidian(PyObject *self, PyObject *args) { substrates_count, stoichiometry, rates, + forms, reactants_lengths, reactants_indexes, From df877d3af1feade7aa6c5d91ad73cc22f855523f Mon Sep 17 00:00:00 2001 From: Heejo Choi Date: Thu, 7 Feb 2019 11:19:54 -0800 Subject: [PATCH 2/7] Add test for optional 'forms' input argument. --- arrow/test/test_arrow.py | 11 +++++++++++ 1 file changed, 11 insertions(+) diff --git a/arrow/test/test_arrow.py b/arrow/test/test_arrow.py index b2f9b83..7ef4516 100644 --- a/arrow/test/test_arrow.py +++ b/arrow/test/test_arrow.py @@ -173,6 +173,17 @@ 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, -1]], np.int64) + + rates = np.array([1, 2]) + forms = np.array([1, 1]) + + arrow = StochasticSystem(stoichiometric_matrix, rates, forms) + result = arrow.evolve(1.0, np.array([50, 20], np.int64)) + if __name__ == '__main__': parser = argparse.ArgumentParser() parser.add_argument('--plot', action='store_true') From 5a60cd9b141e8f3127979b2bdefa9ccc2ccf841f Mon Sep 17 00:00:00 2001 From: Heejo Choi Date: Thu, 7 Feb 2019 16:04:49 -0800 Subject: [PATCH 3/7] Switch propensity computation based on reaction form. --- arrow/obsidian.c | 55 ++++++++++++++++++++++++++++++++-------- arrow/test/test_arrow.py | 9 ++++++- 2 files changed, 52 insertions(+), 12 deletions(-) diff --git a/arrow/obsidian.c b/arrow/obsidian.c index 880e50d..08b17a9 100644 --- a/arrow/obsidian.c +++ b/arrow/obsidian.c @@ -24,6 +24,12 @@ choose(long n, long k) { return combinations; } +// Compute the substrate saturation fraction +double +substrate_fraction(long s, long k) { + return ((double) s / (s + k)); +} + // Perform the Gillespie algorithm with the given stoichiometry, reaction rates and // initial state for the provided duration. @@ -159,18 +165,45 @@ 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[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 with 3 substrates + // Initialize the propensity for reaction with its rate constant + // Note: the rate constant for form '2' reactions is the production of + // k_cat (rate of catalysis) and enzyme count. + propensities[reaction] = rates[reaction]; + + // 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] *= substrate_fraction(count, rates[reaction]); + } + + break; + + default: + printf("arrow.obsidian.evolve - unexpected form: %d", forms[reaction]); } } diff --git a/arrow/test/test_arrow.py b/arrow/test/test_arrow.py index 7ef4516..93259da 100644 --- a/arrow/test/test_arrow.py +++ b/arrow/test/test_arrow.py @@ -182,7 +182,14 @@ def test_forms(): forms = np.array([1, 1]) arrow = StochasticSystem(stoichiometric_matrix, rates, forms) - result = arrow.evolve(1.0, np.array([50, 20], np.int64)) + result = arrow.evolve(10.0, np.array([50, 20], np.int64)) + + print('steps: {}'.format(result['steps'])) + print('time: {}'.format(result['time'])) + print('events: {}'.format(result['events'])) + print('occurrences: {}'.format(result['occurrences'])) + print('outcome: {}'.format(result['outcome'])) + if __name__ == '__main__': parser = argparse.ArgumentParser() From 62d037a16e73f75ae32fd7829410d1ec18a9ed9a Mon Sep 17 00:00:00 2001 From: heejochoi Date: Fri, 8 Feb 2019 13:47:13 -0800 Subject: [PATCH 4/7] Update docstrings. --- arrow/arrow.py | 7 ++++--- arrow/obsidian.c | 30 ++++++++++++++++-------------- 2 files changed, 20 insertions(+), 17 deletions(-) diff --git a/arrow/arrow.py b/arrow/arrow.py index 5b9981c..ef0aba5 100644 --- a/arrow/arrow.py +++ b/arrow/arrow.py @@ -83,9 +83,10 @@ class StochasticSystem(object): def __init__(self, stoichiometry, rates, forms=None, random_seed=0): ''' 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 diff --git a/arrow/obsidian.c b/arrow/obsidian.c index 08b17a9..ded3258 100644 --- a/arrow/obsidian.c +++ b/arrow/obsidian.c @@ -24,9 +24,9 @@ choose(long n, long k) { return combinations; } -// Compute the substrate saturation fraction +// Compute the fractional saturation double -substrate_fraction(long s, long k) { +fraction_saturation(long s, long k) { return ((double) s / (s + k)); } @@ -43,8 +43,10 @@ substrate_fraction(long s, long k) { // * 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. +// * 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. @@ -170,13 +172,15 @@ evolve(MTState *random_state, // 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[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 + // 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]]; @@ -185,19 +189,17 @@ evolve(MTState *random_state, break; - case 1: // Michaelis-Menten with 3 substrates - // Initialize the propensity for reaction with its rate constant - // Note: the rate constant for form '2' reactions is the production of - // k_cat (rate of catalysis) and enzyme count. + case 1: // Michaelis-Menten + // Initialize the propensity for reaction with its maximal propensity propensities[reaction] = rates[reaction]; - // 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. + // 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] *= substrate_fraction(count, rates[reaction]); + propensities[reaction] *= fraction_saturation(count, rates[reaction]); } break; From de698ef9f385c2ea7ec1fa7e9a4340615037a39b Mon Sep 17 00:00:00 2001 From: heejochoi Date: Fri, 8 Feb 2019 21:10:25 -0800 Subject: [PATCH 5/7] Restructure 'rates' into a list of variable length lists. --- arrow/arrow.py | 10 +++++++--- arrow/obsidian.c | 21 ++++++++++++-------- arrow/obsidian.h | 4 +++- arrow/obsidianmodule.c | 44 +++++++++++++++++++++++++++++++----------- 4 files changed, 56 insertions(+), 23 deletions(-) diff --git a/arrow/arrow.py b/arrow/arrow.py index ef0aba5..8ecae92 100644 --- a/arrow/arrow.py +++ b/arrow/arrow.py @@ -102,7 +102,9 @@ def __init__(self, stoichiometry, rates, forms=None, 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) if forms is not None: self.forms = forms @@ -123,7 +125,9 @@ def __init__(self, stoichiometry, rates, forms=None, 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, self.reactants_lengths, self.reactants_indexes, @@ -138,7 +142,7 @@ def __init__(self, stoichiometry, rates, forms=None, 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]) for event in events: occurrences[event] += 1 diff --git a/arrow/obsidian.c b/arrow/obsidian.c index ded3258..dc290d6 100644 --- a/arrow/obsidian.c +++ b/arrow/obsidian.c @@ -40,9 +40,10 @@ fraction_saturation(long s, 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. @@ -91,7 +92,9 @@ 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, @@ -175,7 +178,7 @@ evolve(MTState *random_state, case 0: // Standard Gillespie // Initialize the propensity for reaction with its rate constant - propensities[reaction] = rates[reaction]; + propensities[reaction] = rates_flat[reaction]; // Go through each reactant and calculate its contribution to the // propensity based on the counts of the corresponding substrate, @@ -191,7 +194,8 @@ evolve(MTState *random_state, case 1: // Michaelis-Menten // Initialize the propensity for reaction with its maximal propensity - propensities[reaction] = rates[reaction]; + 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 @@ -199,13 +203,14 @@ evolve(MTState *random_state, for (reactant = 0; reactant < reactants_lengths[reaction]; reactant++) { index = reactants_indexes[reaction] + reactant; count = outcome[reactants[index]]; - propensities[reaction] *= fraction_saturation(count, rates[reaction]); + propensities[reaction] *= fraction_saturation(count, + rates_flat[rates_indexes[reaction] + reactant + 1]); } break; default: - printf("arrow.obsidian.evolve - unexpected form: %d", forms[reaction]); + printf("arrow.obsidian.evolve - unexpected form: %ld", forms[reaction]); } } diff --git a/arrow/obsidian.h b/arrow/obsidian.h index f2f0ea4..ab06354 100644 --- a/arrow/obsidian.h +++ b/arrow/obsidian.h @@ -19,7 +19,9 @@ 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, diff --git a/arrow/obsidianmodule.c b/arrow/obsidianmodule.c index c20752c..2f310fa 100644 --- a/arrow/obsidianmodule.c +++ b/arrow/obsidianmodule.c @@ -39,7 +39,9 @@ typedef struct { int reactions_count; int substrates_count; long *stoichiometry; - double *rates; + double *rates_flat; + long *rates_lengths; + long *rates_indexes; long *forms; long *reactants_lengths; @@ -67,7 +69,9 @@ 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, @@ -105,7 +109,9 @@ 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; @@ -141,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; @@ -190,7 +196,9 @@ 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, @@ -368,7 +376,9 @@ _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, @@ -385,10 +395,12 @@ _invoke_obsidian(PyObject *self, PyObject *args) { *substrates_obj; if (!PyArg_ParseTuple(args, - "iOOOOOOOOOOOOO", + "iOOOOOOOOOOOOOOO", &random_seed, &stoichiometry_obj, - &rates_obj, + &rates_flat_obj, + &rates_lengths_obj, + &rates_indexes_obj, &forms_obj, &reactants_lengths_obj, @@ -412,8 +424,12 @@ _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); @@ -448,7 +464,9 @@ _invoke_obsidian(PyObject *self, PyObject *args) { reactions_count, substrates_count, stoichiometry, - rates, + rates_flat, + rates_lengths, + rates_indexes, forms, reactants_lengths, @@ -470,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); From cf462cce1de667a8d380201fa90a737910c37763 Mon Sep 17 00:00:00 2001 From: heejochoi Date: Fri, 8 Feb 2019 21:14:15 -0800 Subject: [PATCH 6/7] Change constant in fraction saturation computation to double. --- arrow/obsidian.c | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/arrow/obsidian.c b/arrow/obsidian.c index dc290d6..0270bfd 100644 --- a/arrow/obsidian.c +++ b/arrow/obsidian.c @@ -26,7 +26,7 @@ choose(long n, long k) { // Compute the fractional saturation double -fraction_saturation(long s, long k) { +fraction_saturation(long s, double k) { return ((double) s / (s + k)); } From cc2c522fc54aa9cea53a0e2a4033cbcb05d84804 Mon Sep 17 00:00:00 2001 From: heejochoi Date: Fri, 8 Feb 2019 21:15:00 -0800 Subject: [PATCH 7/7] Update test_forms. --- arrow/test/test_arrow.py | 11 ++++++----- 1 file changed, 6 insertions(+), 5 deletions(-) diff --git a/arrow/test/test_arrow.py b/arrow/test/test_arrow.py index 93259da..6f1f130 100644 --- a/arrow/test/test_arrow.py +++ b/arrow/test/test_arrow.py @@ -175,17 +175,18 @@ def test_compare_runtime(): def test_forms(): stoichiometric_matrix = np.array([ - [-1, 1], - [1, -1]], np.int64) + [-1, -1, 1, 0, 0, 0], + [0, 0, 0, -1, -1, 1]], np.int64) - rates = np.array([1, 2]) + 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(10.0, np.array([50, 20], np.int64)) + result = arrow.evolve(1., np.array([50, 50, 0, 20, 20, 0], np.int64)) print('steps: {}'.format(result['steps'])) - print('time: {}'.format(result['time'])) print('events: {}'.format(result['events'])) print('occurrences: {}'.format(result['occurrences'])) print('outcome: {}'.format(result['outcome']))