#294 spatiocyte supports infinite rates

ecell · Sep 27, 2018 · 919f0c5 · 919f0c5
1 parent deaddba
commit 919f0c5
Show file tree

Hide file tree

Showing 3 changed files with 17 additions and 7 deletions.
diff --git a/ecell4/spatiocyte/ReactionEvent.cpp b/ecell4/spatiocyte/ReactionEvent.cpp
@@ -113,14 +113,21 @@ Real FirstOrderReactionEvent::draw_dt()
     const Species& reactant(*(rule_.reactants().begin()));
     const Integer num_r(world_->num_voxels_exact(reactant));
     const Real k(rule_.k());
-    const Real p = k * num_r;
-    Real dt(inf);
-    if (p != 0.)
+    if (num_r > 0)
     {
-        const Real rnd(world_->rng()->uniform(0.,1.));
-        dt = - log(1 - rnd) / p;
+        const Real p = k * num_r;
+        Real dt(inf);
+        if (p != 0.)
+        {
+            const Real rnd(world_->rng()->uniform(0., 1.));
+            dt = - log(1 - rnd) / p;
+        }
+        return dt;
+    }
+    else
+    {
+        return inf;
     }
-    return dt;
 }
 
 } // spatiocyte

diff --git a/ecell4/spatiocyte/StepEvent.cpp b/ecell4/spatiocyte/StepEvent.cpp
@@ -179,7 +179,7 @@ StepEvent::attempt_reaction_(
         const Real k((*itr).k());
         const Real P(k * factor * alpha);
         accp += P;
-        if (accp > 1)
+        if (accp > 1 && !std::isinf(k))
         {
             std::cerr << "The total acceptance probability [" << accp
                 << "] exceeds 1 for '" << speciesA.serial()

diff --git a/ecell4/spatiocyte/utils.cpp b/ecell4/spatiocyte/utils.cpp
@@ -1,3 +1,4 @@
+#include <cmath>
 #include "utils.hpp"
 
 
@@ -69,6 +70,8 @@ const Real calculate_alpha(const ReactionRule& rr, const boost::shared_ptr<Spati
     const ReactionRule::reactant_container_type& reactants(rr.reactants());
     if (reactants.size() != 2)
         return 1.0;
+    else if (std::isinf(rr.k()))
+        return 1.0;
 
     const Species species[2] = {reactants.at(0), reactants.at(1)};
     const MoleculeInfo info[2] = {