Species State Rounding Fix

gillespy2/solvers/cpp/c_base/Tau/tau.cpp

@@ -109,15 +109,15 @@ namespace Gillespy
 		return tau_args;
 	}
 
-	template<typename PType>
+	template<typename PType, typename SType>
 	double select(
 		Gillespy::Model<PType> &model,
 		TauArgs<PType> &tau_args,
 		const double &tau_tol,
 		const double &current_time,
 		const double &save_time,
 		const std::vector<double> &propensity_values,
-		const std::vector<int> &current_state)
+		const std::vector<SType> &current_state)
 	{
 
 		bool critical = false;  // system-wide flag, true when any reaction is critical
@@ -261,14 +261,14 @@ namespace Gillespy
 
 	// Explicitly instantiate initialize/select functions for DISCRETE simulations
 	template TauArgs<double> initialize<double>(Model<double> &model, double tau_tol);
-	template double select<double>(
+	template double select<double, double>(
 		Model<double> &model,
 		TauArgs<double> &tau_args,
 		const double &tau_tol,
 		const double &current_time,
 		const double &save_time,
 		const std::vector<double> &propensity_values,
-		const std::vector<int> &current_state);
+		const std::vector<double> &current_state);
 
 	// Explicitly instantiate initialize/select functions for HYBRID simulations
 	template TauArgs<unsigned int> initialize<unsigned int>(Model<unsigned int> &model, double tau_tol);

gillespy2/solvers/cpp/c_base/Tau/tau.h

@@ -47,27 +47,27 @@ namespace Gillespy
     template<typename PType>
     TauArgs<PType> initialize(Model<PType> &model, double tau_tol);
 
-    template<typename PType>
+    template<typename PType, typename SType = int>
     double select(
         Model<PType> &model,
         TauArgs<PType> &tau_args,
         const double &tau_tol,
         const double &current_time,
         const double &save_time,
         const std::vector<double> &propensity_values,
-        const std::vector<int> &current_state);
+        const std::vector<SType> &current_state);
 
     // Continuous tau args is used by hybrid solver
     template struct TauArgs<double>;
     extern template TauArgs<double> initialize(Model<double> &model, double tau_tol);
-    extern template double select<double>(
+    extern template double select<double, double>(
         Model<double> &model,
         TauArgs<double> &tau_args,
         const double &tau_tol,
         const double &current_time,
         const double &save_time,
         const std::vector<double> &propensity_values,
-        const std::vector<int> &current_state);
+        const std::vector<double> &current_state);
 
     // Discrete tau args is used by tau leaping solver
     template struct TauArgs<unsigned int>;

gillespy2/solvers/cpp/c_base/tau_hybrid_cpp_solver/HybridModel.cpp

@@ -203,11 +203,19 @@ namespace Gillespy
 
 		HybridReaction::HybridReaction()
 				: mode(SimulationState::DISCRETE),
-				  base_reaction(nullptr)
+				  m_base_reaction(nullptr),
+				  m_id(-1)
 		{
 			// Empty constructor body
 		}
 
+		HybridReaction::HybridReaction(Reaction *base_reaction)
+			: HybridReaction()
+		{
+			base_reaction = base_reaction;
+			m_id = base_reaction->id;
+		}
+
 		HybridSpecies::HybridSpecies()
 				: user_mode(SimulationState::DYNAMIC),
 				  partition_mode(SimulationState::DISCRETE),
@@ -235,15 +243,14 @@ namespace Gillespy
 
 			for (int rxn_i = 0; rxn_i < model.number_reactions; ++rxn_i)
 			{
-				reaction_state[rxn_i].base_reaction = &model.reactions[rxn_i];
+				reaction_state[rxn_i].set_base_reaction(&model.reactions[rxn_i]);
 			}
 		}
 
 
 		double DifferentialEquation::evaluate(
 				const double t,
-				double *ode_state,
-				int *ssa_state)
+				double *ode_state)
 		{
 			double sum = 0.0;
 
@@ -254,7 +261,7 @@ namespace Gillespy
 
 			for (auto &formula : formulas)
 			{
-				sum += formula(ode_state, ssa_state);
+				sum += formula(ode_state);
 			}
 
 			return sum;
@@ -273,9 +280,8 @@ namespace Gillespy
 				spec.diff_equation.formulas.clear();
 			}
 
-			for (int rxn_i = 0; rxn_i < reactions.size(); ++rxn_i)
+			for (HybridReaction &rxn : reactions)
 			{
-				HybridReaction rxn = reactions[rxn_i];
 				if (rxn.mode == SimulationState::DISCRETE)
 				{
 					continue;
@@ -284,38 +290,18 @@ namespace Gillespy
 				for (int spec_i = 0; spec_i < species.size(); ++spec_i)
 				{
 					// A species change of 0 indicates that this species is not a dependency for this reaction.
-					if (rxn.base_reaction->species_change[spec_i] == 0)
+					if (rxn.get_base_reaction()->species_change[spec_i] == 0)
 					{
 						continue;
 					}
 
 					HybridSpecies &spec = species[spec_i];
 					auto &formula_set = spec.diff_equation.formulas;
-					int spec_diff = rxn.base_reaction->species_change[spec_i];
+					int spec_diff = rxn.get_base_reaction()->species_change[spec_i];
 
-					switch (spec.partition_mode)
-					{
-					case SimulationState::CONTINUOUS:
-						formula_set.push_back([rxn_i, spec_diff](
-								double *ode_state,
-								int *ssa_state)
-											  {
-												  return spec_diff * Reaction::propensity(rxn_i, ode_state);
-											  });
-						break;
-
-					case SimulationState::DISCRETE:
-						formula_set.push_back([rxn_i, spec_diff](
-								double *ode_state,
-								int *ssa_state)
-											  {
-												  return spec_diff * Reaction::propensity(rxn_i, ssa_state);
-											  });
-						break;
-
-					default:
-						break;
-					}
+					formula_set.emplace_back([&rxn, spec_diff](double *state) {
+						return spec_diff * rxn.propensity(state);
+					});
 				}
 			}
 		}
@@ -342,7 +328,7 @@ namespace Gillespy
 				{
 					// Reaction has a dependency on a species if its dx is positive or negative.
 					// Any species with "dependency" change of 0 is by definition not a dependency.
-					if (rxn.base_reaction->species_change[spec_i] == 0)
+					if (rxn.get_base_reaction()->species_change[spec_i] == 0)
 					{
 						continue;
 					}
@@ -410,7 +396,7 @@ namespace Gillespy
 					// Selected species is either a reactant or a product, depending on whether
 					//   dx is positive or negative.
 					// 0-dx species are not dependencies of this reaction, so dx == 0 is ignored.
-					int spec_dx = rxn.base_reaction->species_change[spec_i];
+					int spec_dx = rxn.get_base_reaction()->species_change[spec_i];
 					if (spec_dx < 0)
 					{
 						// Selected species is a reactant.

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

@@ -177,10 +177,10 @@ namespace Gillespy
 		struct DifferentialEquation
 		{
 		public:
-			std::vector<std::function<double(double *, int *)>> formulas;
+			std::vector<std::function<double(double *)>> formulas;
 			std::vector<std::function<double(double, double *, double *, const double *)>> rate_rules;
 
-			double evaluate(double t, double *ode_state, int *ssa_state);
+			double evaluate(double t, double *ode_state);
 		};
 
 		enum SimulationState : unsigned int
@@ -227,10 +227,62 @@ namespace Gillespy
 
 		struct HybridReaction
 		{
-			Reaction *base_reaction;
 			SimulationState mode;
 
+			inline double propensity(double *state) const
+			{
+				return propensity(state, Reaction::s_variables.get(), Reaction::s_constants.get());
+			}
+
+			inline double propensity(double *state, double *parameters, const double *constants) const
+			{
+				switch (mode)
+				{
+				case SimulationState::CONTINUOUS:
+					return ode_propensity(state, parameters, constants);
+				case SimulationState::DISCRETE:
+				default:
+					return ssa_propensity(state, parameters, constants);
+				}
+			}
+
+			inline double ode_propensity(double *state) const
+			{
+				return ode_propensity(state, Reaction::s_variables.get(), Reaction::s_constants.get());
+			}
+
+			inline double ode_propensity(double *state, double *parameters, const double *constants) const
+			{
+				return map_ode_propensity(m_id, state, parameters, constants);
+			}
+
+			inline double ssa_propensity(double *state) const
+			{
+				return ssa_propensity(state, Reaction::s_variables.get(), Reaction::s_constants.get());
+			}
+
+			inline double ssa_propensity(double *state, double *parameters, const double *constants) const
+			{
+				return map_ssa_propensity(m_id, state, parameters, constants);
+			}
+
+			inline void set_base_reaction(Reaction *reaction)
+			{
+				m_base_reaction = reaction;
+				m_id = reaction->id;
+			}
+
+			inline Reaction *get_base_reaction() const
+			{
+				return m_base_reaction;
+			}
+
 			HybridReaction();
+			explicit HybridReaction(Reaction *base_reaction);
+
+		private:
+			Reaction *m_base_reaction;
+			ReactionId m_id;
 		};
 
 		struct HybridSimulation : Simulation<double>

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

@@ -87,14 +87,12 @@ namespace Gillespy
 				// TODO: change back double -> hybrid_state, once we figure out how that works
 				EventList event_list;
 				std::vector<double> current_state(num_species);
-				std::vector<int> current_populations(num_species);
 
 				// Initialize the species population for the trajectory.
 				for (int spec_i = 0; spec_i < num_species; ++spec_i)
 				{
 					current_state[spec_i] = species[spec_i].initial_population;
 					simulation->current_state[spec_i] = current_state[spec_i];
-					current_populations[spec_i] = species[spec_i].initial_population;
 				}
 				simulation->reset_output_buffer(traj);
 				simulation->output_buffer_range(std::cout);
@@ -143,32 +141,20 @@ namespace Gillespy
 					for (unsigned int rxn_i = 0; rxn_i < num_reactions; ++rxn_i)
 					{
 						HybridReaction &rxn = simulation->reaction_state[rxn_i];
-						double propensity = 0.0;
-						switch (rxn.mode)
-						{
-						case SimulationState::CONTINUOUS:
-							propensity = Reaction::propensity(rxn_i, current_state.data());
-							break;
-						case SimulationState::DISCRETE:
-							propensity = Reaction::propensity(rxn_i, current_populations.data());
-							break;
-						default:
-							break;
-						}
-						sol.data.propensities[rxn_i] = propensity;
+						sol.data.propensities[rxn_i] = rxn.propensity(current_state.data());
 					}
 					if (interrupted)
 						break;
 
 					// Expected tau step is determined.
-					tau_step = select(
+					tau_step = select<double, double>(
 							model,
 							tau_args,
 							tau_tol,
 							simulation->current_time,
 							save_time,
 							sol.data.propensities,
-							current_populations
+							current_state
 					);
 					partition_species(
 							simulation->reaction_state,
@@ -304,7 +290,6 @@ namespace Gillespy
 									current_state[p_i] += population_changes[p_i];
 									result.concentrations[p_i] = current_state[p_i];
 								}
-								current_populations[p_i] = (int) current_state[p_i];
 							}
 						}
 					} while (invalid_state && !interrupted);