Merge pull request #22 from CovertLab/standardization

Language and mathematical standardizations

README.md

@@ -1,62 +1,81 @@
-# arrow
+# Arrow
 
-“... even if the previous millisecond is closer to us than the birth of the universe, it is equally out of reach.” 
+“... even if the previous millisecond is closer to us than the birth of the universe, it is equally out of reach.”
 ― Jean-Christophe Valtat, Luminous Chaos
 
-## concept
+## Concept
 
-This library implements a generalized version of the [Gillespie Algorithm](https://en.wikipedia.org/wiki/Gillespie_algorithm), a stochastic approach to numerically solving discrete systems. Each iteration, the algorithm will calculate the propensities for each reaction given a rate and the counts of the reactants present in the current state of the system, then selects one reaction to occur and the interval of time between the previous reaction and the current reaction. Iterating this produces a trajectory (or `history`) of the state vector over the course of the simulation.
+This library implements a generalized version of the [Gillespie
+Algorithm](https://en.wikipedia.org/wiki/Gillespie_algorithm), a stochastic
+approach to numerically solving discrete systems. Each iteration, the algorithm
+will calculate the propensities for each reaction given a rate and the counts
+of the reactants present in the current state of the system, then selects one
+reaction to occur and the interval of time between the previous reaction and
+the current reaction. Iterating this produces a trajectory (or `history`) of
+the state vector over the course of the simulation.
 
-## installation
+## Installation
 
-Add the following to your `requirements.txt`, or `pip install stochastic-arrow`:
+Add the following to your `requirements.txt`, or
+`pip install stochastic-arrow`:
 
     stochastic-arrow==0.0.1
 
-## usage
+## Usage
 
-The `arrow` library presents a single class as an interface, `StochasticSystem`, which operates on a set of reactions (encoded as a `numpy` matrix) and associated reaction rates:
+The `arrow` library presents a single class as an interface,
+`StochasticSystem`, which operates on a set of reactions (encoded as a `numpy`
+matrix of stoichiometrix coefficients) and associated reaction rates:
 
 ```python
 from arrow import StochasticSystem
 import numpy as np
 
-# Each row is a reaction and each column is a molecular species (or other entity).
-# The top reaction here means that the first and second elements combine to create the third,
-# while the fourth is unaffected.
-reactions = np.array([
-    [-1, -1, 1, 0],
-    [-2, 0, 0, 1],
-    [1, 1, -1, 0]])
+# Each column is a reaction and each row is a molecular species (or other
+# entity). The first reaction here means that the first and second elements
+# combine to create the third, while the fourth is unaffected.
+stoichiometric_matrix = np.array([
+    [-1, -2, +1],
+    [-1,  0, +1],
+    [+1,  0, -1],
+    [ 0, +1,  0]
+    ])
 
 # Each reaction has an associated rate for how probable that reaction is.
-rates = np.array([3, 1, 1])
+rates = np.array([3.0, 1.0, 1.0])
 
-# Once we have a matrix of reactions and their associated rates, we can construct the system.
-system = StochasticSystem(reactions, rates)
+# Once we have a matrix of reactions and their associated rates, we can
+# construct the system.
+system = StochasticSystem(stoichiometric_matrix, rates)
 ```
 
-Now that the system has been instantiated, we can invoke it with any initial state vector and then run it for a given time interval:
+Now that the system has been instantiated, we can invoke it with any initial
+state vector and then run it for a given time interval:
 
 ```python
-# This gives the initial state of the system (counts of each molecular species, for instance)
+# This gives the initial state of the system (counts of each molecular species,
+# for instance).
 state = np.array([1000, 1000, 0, 0])
 
-# We also specify how long we want the simulation to run. Here we set it to one second
+# We also specify how long we want the simulation to run. Here we set it to one
+# second.
 duration = 1
 
-# Once we have an initial state and duration, we can run the simulation for the given duration.
-# `evolve` returns the history of the state vector for each time step, and the history of time
-# steps as they will be in uneven increments throughout the simulation.
-history, steps = system.evolve(state, duration)
+# Once we have an initial state and duration, we can run the simulation for the
+# given duration. `evolve` returns the history of the state vector (counts) for
+# each time step, and the history of time steps as they will be in uneven
+# increments throughout the simulation.
+time, counts = system.evolve(state, duration)
 ```
 
-## testing
+## Testing
 
-`arrow` uses pytest: https://docs.pytest.org/en/latest/ so you can test simply by saying:
+`arrow` uses pytest: https://docs.pytest.org/en/latest/ so you can test simply
+by invoking:
 
     > pytest
 
-Also, we have a test that generates plots of various systems which can be run like so:
+Also, we have a test that generates plots of various systems which can be run
+like so:
 
-    > python arrow/test/test_arrow.py
+    > python arrow/test/test_arrow.py

arrow/arrow.py

@@ -9,86 +9,93 @@ def choose(n, k):
     return np.rint(product)
 
 
-def propensity(reaction, state, form):
-    reactants = np.where(reaction < 0)
+def propensity(stoichiometric_matrix, state, form):
+    reactants = np.where(stoichiometric_matrix < 0)
     terms = [
-        form(state[reactant], -reaction[reactant])
+        form(state[reactant], -stoichiometric_matrix[reactant])
         for reactant in reactants[0]]
 
     return np.array(terms).prod()
 
-
-def step(reactions, rates, state, forms, propensities=[], update_reactions=()):
+def step(stoichiometric_matrix, rates, state, forms, propensities=[], update_reactions=()):
     if len(update_reactions):
         for update in update_reactions:
-            reaction = reactions[update]
+            stoichiometry = stoichiometric_matrix[:, update]
             form = forms if callable(forms) else forms[update]
-            propensities[update] = propensity(reaction, state, form)
+            propensities[update] = propensity(stoichiometry, state, form)
     else:
-        propensities = [
-            propensity(reaction, state, forms if callable(forms) else forms[index])
-            for index, reaction in enumerate(reactions)]
+        propensities = np.array([
+            propensity(stoichiometry, state, forms if callable(forms) else forms[index])
+            for index, stoichiometry in enumerate(stoichiometric_matrix.T)])
 
     distribution = (rates * propensities)
     total = distribution.sum()
-    if not total:
-        return (state, 0, -1, propensities)
-
-    dt = np.random.exponential(1 / total)
-    random = np.random.uniform(0, 1) * total
-
-    progress = 0
-    for choice, interval in enumerate(distribution):
-        progress += interval
-        if random <= progress:
-            break
 
-    reaction = reactions[choice]
-    outcome = state + reaction
+    if total == 0:
+        time_to_next = 0
+        outcome = state
+        choice = -1
+
+    else:
+        time_to_next = np.random.exponential(1 / total)
+        random = np.random.uniform(0, 1) * total
+
+        progress = 0
+        for choice, interval in enumerate(distribution):
+            progress += interval
+            if random <= progress:
+                break
+
+        stoichiometry = stoichiometric_matrix[:, choice]
+        outcome = state + stoichiometry
 
-    return outcome, dt, choice, propensities
+    return time_to_next, outcome, choice, propensities
 
 
-def evolve(reactions, rates, state, duration, forms=choose):
-    time = 0
-    history = [state]
-    steps = [0]
+def evolve(stoichiometric_matrix, rates, state, duration, forms=choose):
+    time_current = 0
+    time = [0]
+    counts = [state]
     propensities = []
     update_reactions = []
     events = np.zeros(rates.shape)
 
     while True:
-        state, dt, choice, propensities = step(
-            reactions,
+        time_to_next, state, choice, propensities = step(
+            stoichiometric_matrix,
             rates,
             state,
             forms,
             propensities,
             update_reactions)
 
-        time += dt
-        if not dt or time > duration:
+        time_current += time_to_next
+        if not time_to_next or time_current > duration:
             break
 
-        history.append(state)
-        steps.append(time)
+        time.append(time_current)
+        counts.append(state)
 
         events[choice] += 1
-        reaction = reactions[choice]
-        involved = np.where(reaction != 0)
-        update_reactions = np.where(reactions[:, involved] != 0)[0]
 
-    return history, steps, events
+        stoichiometry = stoichiometric_matrix[:, choice]
+        involved = np.where(stoichiometry != 0)[0]
+        update_reactions = np.where(np.any(stoichiometric_matrix[involved] < 0, 0))[0]
+
+    time = np.array(time)
+    counts = np.array(counts)
+
+    return time, counts, events
 
 
 class StochasticSystem(object):
-    def __init__(self, reactions, rates, forms=None):
-        self.reactions = reactions
+    def __init__(self, stoichiometric_matrix, rates, forms=None):
+        self.stoichiometric_matrix = stoichiometric_matrix
         self.rates = rates
         self.forms = forms or choose
 
     def step(self, state):
-        return step(self.reactions, self.rates, state, forms=self.forms)
+        return step(self.stoichiometric_matrix, self.rates, state, forms=self.forms)
 
     def evolve(self, state, duration):
-        return evolve(self.reactions, self.rates, state, duration, forms=self.forms)
+        return evolve(self.stoichiometric_matrix, self.rates, state, duration, forms=self.forms)