Merge pull request #860 from StochSS/develop

Release v1.7.1

AUTHORS

@@ -4,7 +4,8 @@ Ghilman Brock
 Eliot Dixon
 Dalton Nickerson
 Sean Matthew
-George Hall
+Matthew Geiger
+Kassidy Hall
 W.R. Jackson
 Samuel Hodges
 Emma Weisgerber
@@ -16,3 +17,7 @@ Bryan Rumsey
 Mason Kidwell
 Jesse Reeve
 Fin Carter
+Ethan Green
+Josh Cooper
+Matthew Dippel
+Joshua Fisher-Caudill

README.md

@@ -176,6 +176,7 @@ Authors and history
 * [**Eliot Dixon**](https://github.com/edixon1)
 * [**Dalton Nickerson**](https://github.com/Thisisnotdalton)
 * [**Sean Matthew**](https://github.com/seanebum)
+* [**Matthew Geiger**](https://github.com/popensesame)
 * [**Kassidy Hall** ](https://github.com/technodivergent)
 * [**W.R. Jackson** ](https://github.com/JustJackson)
 * [**Samuel Hodges**](https://github.com/hodgespodge)
@@ -188,6 +189,11 @@ Authors and history
 * [**Mason Kidwell**](https://github.com/makdl)
 * [**Jesse Reeve**](https://github.com/jdreeve)
 * [**Fin Carter**](https://github.com/Fin109)
+* [**Ethan Green**](https://github.com/ethangreen-dev)
+* [**Josh Cooper**](https://github.com/jtcooper10)
+* [**Matthew Dippel**](https://github.com/mdip226)
+* [**Joshua Fisher-Caudill**](https://github.com/joshuafisherc)
+
 
 Acknowledgments
 ---------------

gillespy2/__version__.py

@@ -22,7 +22,7 @@
 # =============================================================================
 
 
-__version__      = '1.7.0'
+__version__      = '1.7.1'
 
 __title__        = 'GillesPy2'
 __description__  = 'Python interface for Gillespie-style biochemical simulations'
@@ -31,4 +31,4 @@
 __author__       = 'See AUTHORS'
 __email__        = '[email protected]'
 __license__      = 'GPL'
-__copyright__    = 'Copyright (C) 2017-2021'
+__copyright__    = 'Copyright (C) 2017-2022'

gillespy2/core/events.py

@@ -39,12 +39,15 @@ class EventAssignment(Jsonify):
 
     """
 
+
     def __init__(self, name=None, variable=None, expression=None):
 
         if name in (None, ""):
-            self.name = f'evn{uuid.uuid4()}'.replace('-', '_')
+            name = f'evn{uuid.uuid4()}'.replace('-', '_')
         else:
-            self.name = name
+            from gillespy2.core import log
+            log.warning("EventAssignment.name has been deprecated.")
+        self.__name_deprecated = name
 
         self.variable = variable
         self.expression = expression
@@ -70,6 +73,14 @@ def __init__(self, name=None, variable=None, expression=None):
     def __str__(self):
         return f"{self.variable}: {self.expression}"
 
+    def __getattr__(self, key):
+        if key == 'name':
+            from gillespy2.core import log
+            log.warning('EventAssignment.name has been deprecated.')
+            return self.__name_deprecated
+        else: 
+            raise AttributeError
+
 class EventTrigger(Jsonify):
     """
     Trigger detects changes in model/environment conditions in order to fire an

gillespy2/core/model.py

@@ -1062,6 +1062,9 @@ def make_translation_table(self):
         assignments = self.listOfAssignmentRules.values()
         rates = self.listOfRateRules.values()
         events = self.listOfEvents.values()
+        for event in events:
+            for assignment in event.__dict__['assignments']:
+                del assignment.__dict__['__name_deprecated']
         functions = self.listOfFunctionDefinitions.values()
 
         # A translation table is used to anonymize user-defined variable names and formulas into generic counterparts.

gillespy2/core/reaction.py

@@ -295,7 +295,7 @@ def _create_mass_action(self):
                 reactant = reactant.name
             # Case 1: 2X -> Y
             if stoichiometry == 2:
-                propensity_function = f"0.5 * {propensity_function} * {reactant} * ({reactant} - 1) / vol"
+                propensity_function = f"{propensity_function} * {reactant} * ({reactant} - 1) / vol"
                 ode_propensity_function += f" * {reactant} * {reactant}"
             else:
                 # Case 3: X1, X2 -> Y;

gillespy2/solvers/cpp/c_base/arg_parser.cpp

@@ -95,10 +95,12 @@ char ArgParser::match_arg(std::string &token)
         return 'M';
     }
 
-    else
+    if (!token.compare("--use_root_finding"))
     {
-        return 0;
+        return 'u';
     }
+
+    return 0;
 }
 
 ArgParser::ArgParser(int argc, char *argv[])
@@ -176,6 +178,10 @@ ArgParser::ArgParser(int argc, char *argv[])
                 verbose = true;
                 break;
 
+            case 'u':
+                use_root_finding = true;
+                break;
+
             case 'R':
                 std::stringstream(argv[i + 1]) >> rtol;
                 break;

gillespy2/solvers/cpp/c_base/arg_parser.h

@@ -51,6 +51,7 @@ class ArgParser
     double atol = 1e-12;
 
     bool verbose = false;
+    bool use_root_finding = false;
 
     ArgParser(int argc, char *argv[]);
     ~ArgParser();

gillespy2/solvers/cpp/c_base/tau_hybrid_cpp_solver/TauHybridSimulation.cpp

@@ -78,9 +78,10 @@ int main(int argc, char* argv[])
         parser.rtol,
         parser.atol,
         parser.max_step,
+
     };
 
-    TauHybrid::TauHybridCSolver(&simulation, events, logger, tau_tol, config);
+    TauHybrid::TauHybridCSolver(&simulation, events, logger, tau_tol, config, parser.use_root_finding);
     simulation.output_buffer_final(std::cout);
     return simulation.get_status();
 }

gillespy2/solvers/cpp/c_base/tau_hybrid_cpp_solver/TauHybridSolver.cpp

@@ -84,6 +84,7 @@ namespace Gillespy
                     // Temporary variable for the reaction's state.
                     // Does not get updated unless the changes are deemed valid.
                     double rxn_state = result.reactions[rxn_i];
+                    double old_rxn_state = rxn_state;
 
                     if (simulation->reaction_state[rxn_i].mode == SimulationState::DISCRETE) {
                         unsigned int rxn_count = 0;
@@ -109,15 +110,14 @@ namespace Gillespy
             }
         }
 
-        bool TakeIntegrationStep(Integrator&sol, IntegrationResults&result, double next_time, int*population_changes,
+        bool TakeIntegrationStep(Integrator&sol, IntegrationResults&result, double *next_time, int*population_changes,
          std::vector<double> current_state, std::set<unsigned int>&rxn_roots, 
          std::set<int>&event_roots, HybridSimulation*simulation, URNGenerator&urn, 
          int only_reaction_to_fire){
             // Integration Step
             // For deterministic reactions, the concentrations are updated directly.
             // For stochastic reactions, integration updates the rxn_offsets vector.
-            //IntegrationResults result = sol.integrate(&next_time, event_roots, rxn_roots);
-            result = sol.integrate(&next_time, event_roots, rxn_roots);
+            result = sol.integrate(next_time, event_roots, rxn_roots);
             if (sol.status == IntegrationStatus::BAD_STEP_SIZE)
             {
                 simulation->set_status(HybridSimulation::INTEGRATOR_FAILED);
@@ -134,10 +134,11 @@ namespace Gillespy
 
 
 
-        bool IsStateNegativeCheck(int num_species, int*population_changes, std::vector<double> current_state){
+        bool IsStateNegativeCheck(int num_species, int*population_changes, std::vector<double> current_state, std::set<Species<double>>tau_args_reactants){  
             // Explicitly check for invalid population state, now that changes have been tallied.
-            for (int spec_i = 0; spec_i < num_species; ++spec_i) {
-                if (current_state[spec_i] + population_changes[spec_i] < 0) {
+            // Note: this should only check species that are reactants or products
+            for (const auto &r : tau_args_reactants) {
+                if (population_changes[r.id] != 0 && current_state[r.id] + population_changes[r.id] < 0) {
                     return true;
                 }
             }
@@ -150,7 +151,8 @@ namespace Gillespy
                 std::vector<Event> &events,
                 Logger &logger,
                 double tau_tol,
-                SolverConfiguration config)
+                SolverConfiguration config,
+                bool default_use_root_finding)
         {
             if (simulation == NULL)
             {
@@ -163,6 +165,9 @@ namespace Gillespy
             int num_reactions = model.number_reactions;
             int num_trajectories = simulation->number_trajectories;
             std::unique_ptr<Species<double>[]> &species = model.species;
+            bool use_root_finding = default_use_root_finding;
+            bool in_event_handling = false;
+            unsigned int neg_state_loop_cnt = 0;
 
             generator = std::mt19937_64(simulation->random_seed);
             URNGenerator urn(simulation->random_seed);
@@ -319,101 +324,58 @@ namespace Gillespy
 
                     IntegrationResults result;
 
-                    if(!TauHybrid::TakeIntegrationStep(sol, result, next_time, population_changes, current_state, rxn_roots, event_roots, simulation, urn, -1)){
+                    if(in_event_handling){
+                        sol.use_events(events, simulation->reaction_state);
+                        sol.enable_root_finder();
+                    }else if(use_root_finding){
+                        sol.use_reactions(simulation->reaction_state);
+                        sol.enable_root_finder();
+                        if(neg_state_loop_cnt > 0){
+                            neg_state_loop_cnt--;
+                        }else{
+                            use_root_finding = default_use_root_finding;
+                        }
+                    }else{
+                        sol.disable_root_finder();
+                    }
+
+
+                    if(!TauHybrid::TakeIntegrationStep(sol, result, &next_time, population_changes, current_state, rxn_roots, event_roots, simulation, urn, -1)){
                         return;
                     }
 
                     // Check if we have gone negative
-                    if (TauHybrid::IsStateNegativeCheck(num_species, population_changes, current_state)) {
+                    if (TauHybrid::IsStateNegativeCheck(num_species, population_changes, current_state, tau_args.reactants)) {
                         // If state is invalid, we took too agressive tau step and need to take a single SSA step forward
                         // Restore the solver to the intial step state
                         sol.restore_state();
-
-                        // Calculate floor()'ed state for use in SSA step
-                        for(int spec_i = 0; spec_i < num_species; ++spec_i){
-                            floored_current_state[spec_i] = floor(current_state[spec_i]);
-                        }
-                        // estimate the time to the first stochatic reaction by assuming constant propensities
-                        double min_tau = 0.0;
-                        int rxn_selected = -1;
-
-                        double *rxn_state = sol.get_reaction_state();
-
-                        for (int rxn_k = 0; rxn_k < num_reactions; ++rxn_k) {
-                            HybridReaction &rxn = simulation->reaction_state[rxn_k];
-                            double propensity_value = rxn.ssa_propensity(current_state.data());
-                            double floored_propensity_value = rxn.ssa_propensity(floored_current_state);
-                            //estimate the zero crossing time
-                            if(floored_propensity_value > 0.0){
-                                double est_tau = -1*  rxn_state[rxn_k] / propensity_value;
-
-                                if(rxn_selected == -1 || est_tau < min_tau ){
-                                    min_tau = est_tau;
-                                    rxn_selected = rxn_k;
-                                }
-                            }
-                        }
-                        if(rxn_selected == -1){
-                            simulation->set_status(HybridSimulation::NEGATIVE_STATE_NO_SSA_REACTION);
-                            return;
-                        }
-                        // if min_tau < 1e-10, we can't take an ODE step that small.
-                        if( min_tau < 1e-10 ){
-                            // instead we will fire the reaction
-                            CalculateSpeciesChangeAfterStep(result, population_changes, current_state, rxn_roots,  event_roots, simulation, urn, rxn_selected);
-                            // re-attempt the step at the same time 
-                            next_time = simulation->current_time;
-
-                        }else{
-                            // Use the found tau-step for single SSA
-                            next_time = simulation->current_time + min_tau;
-
-                            //***********************************
-                            //***********************************
-
-                            // Integreate the system forward
-                            if(!TauHybrid::TakeIntegrationStep(sol, result, next_time, population_changes, current_state, rxn_roots,  event_roots, simulation, urn, rxn_selected)){
-                                return;
-                            }
-                        }
-                        // check for invalid state again
-                        if (TauHybrid::IsStateNegativeCheck(num_species, population_changes, current_state)) {
-                            //Got an invalid state after the SSA step
-                            simulation->set_status(HybridSimulation::INVALID_AFTER_SSA);
-                            return;
-                        }
+                        use_root_finding=true;
+                        neg_state_loop_cnt = 2; // How many single SSA events should we find before we go back to tau steping
+                        continue;
                     }
 
-                    // "Permanently" update the rxn_state and populations.
+
+                    // Update solver object with stochastic changes
                     for (int p_i = 0; p_i < num_species; ++p_i)
                     {
                         if (!simulation->species_state[p_i].boundary_condition)
                         {
-                            // Boundary conditions are not modified directly by reactions.
-                            // As such, population dx in stochastic regime is not considered.
-                            // For deterministic species, their effective dy/dt should always be 0.
                             HybridSpecies *spec = &simulation->species_state[p_i];
                             if( spec->partition_mode == SimulationState::CONTINUOUS ){
-                                current_state[p_i] = result.concentrations[p_i] + population_changes[p_i];
+                                result.concentrations[p_i] = result.concentrations[p_i] + population_changes[p_i];
                             }else if( spec->partition_mode == SimulationState::DISCRETE ){
-                                current_state[p_i] += population_changes[p_i];
+                                result.concentrations[p_i] = current_state[p_i] + population_changes[p_i];
                             }
-                            result.concentrations[p_i] = current_state[p_i]; 
                         }
                     }
-
-                    if (interrupted){
-                        break;
-                    }
-
                     // ===== <EVENT HANDLING> =====
                     if (!event_list.has_active_events())
                     {
                         if (event_list.evaluate_triggers(N_VGetArrayPointer(sol.y), next_time))
                         {
                             sol.restore_state();
-                            sol.use_events(events, simulation->reaction_state);
-                            sol.enable_root_finder();
+                            use_root_finding=true;
+                            in_event_handling=true;
                             continue;
                         }
                     }
@@ -422,12 +384,31 @@ namespace Gillespy
                         double *event_state = N_VGetArrayPointer(sol.y);
                         if (!event_list.evaluate(event_state, num_species, next_time, event_roots))
                         {
-                            sol.disable_root_finder();
+                            in_event_handling=false;
+                            use_root_finding = default_use_root_finding; // set to default
                         }
                         std::copy(event_state, event_state + num_species, current_state.begin());
                     }
                     // ===== </EVENT HANDLING> =====
 
+                    // "Permanently" update species populations.
+                    // (needs to be below event handling)
+                    for (int p_i = 0; p_i < num_species; ++p_i)
+                    {
+                        if (!simulation->species_state[p_i].boundary_condition)
+                        {
+                            // Boundary conditions are not modified directly by reactions.
+                            // As such, population dx in stochastic regime is not considered.
+                            // For deterministic species, their effective dy/dt should always be 0.
+                            current_state[p_i]  = result.concentrations[p_i];
+                        }
+                    }
+
+                    if (interrupted){
+                        break;
+                    }
+
+
                     // Output the results for this time step.
                     sol.refresh_state();
                     simulation->current_time = next_time;

gillespy2/solvers/cpp/c_base/tau_hybrid_cpp_solver/TauHybridSolver.h

@@ -25,6 +25,6 @@ namespace Gillespy
 {
     namespace TauHybrid
     {
-        void TauHybridCSolver(HybridSimulation* simulation, std::vector<Event> &events, Logger &logger, double tau_tol, SolverConfiguration config);
+        void TauHybridCSolver(HybridSimulation* simulation, std::vector<Event> &events, Logger &logger, double tau_tol, SolverConfiguration config, bool use_root_finding);
     }
 }