🐛 Fix reactions with OffLatticeSpace

ecell · Jan 23, 2020 · 0f3e1c4 · 0f3e1c4
1 parent b0fe1aa
commit 0f3e1c4
Show file tree

Hide file tree

Showing 4 changed files with 27 additions and 21 deletions.
diff --git a/ecell4/spatiocyte/SpatiocyteReactions.cpp b/ecell4/spatiocyte/SpatiocyteReactions.cpp
@@ -12,8 +12,7 @@ namespace spatiocyte
 inline const std::string get_serial(boost::shared_ptr<SpatiocyteWorld> world,
                                     const Voxel &voxel)
 {
-    boost::shared_ptr<const VoxelPool> mtype(voxel.get_voxel_pool());
-    return mtype->is_vacant() ? "" : mtype->species().serial();
+    return voxel.get_voxel_pool()->species().serial();
 }
 
 inline const std::string get_location(boost::shared_ptr<SpatiocyteWorld> world,
@@ -22,8 +21,7 @@ inline const std::string get_location(boost::shared_ptr<SpatiocyteWorld> world,
     boost::shared_ptr<const VoxelPool> mtype(voxel.get_voxel_pool());
     if (mtype->is_vacant())
         return "";
-    boost::shared_ptr<const VoxelPool> ltype(mtype->location());
-    return ltype->species().serial();
+    return mtype->location()->species().serial();
 }
 
 static inline void make_product(boost::shared_ptr<SpatiocyteWorld> world,
@@ -444,7 +442,7 @@ apply_second_order_reaction(boost::shared_ptr<SpatiocyteWorld> world,
         return apply_vanishment(world, reactant_item0, reactant_item1);
     case 1:
         return apply_ab2c(world, reactant_item0, reactant_item1,
-                          *(products.begin()));
+                          products.at(0));
     case 2:
         return apply_ab2cd(world, reactant_item0, reactant_item1,
                            products.at(0), products.at(1));

diff --git a/ecell4/spatiocyte/SpatiocyteWorld.cpp b/ecell4/spatiocyte/SpatiocyteWorld.cpp
@@ -361,6 +361,17 @@ boost::optional<Voxel> SpatiocyteWorld::check_neighbor(const Voxel &voxel,
         }
     }
 
+    if (const auto neighbors = neighbors_.find(voxel))
+    {
+        for (const auto &neighbor : *neighbors)
+        {
+            if (neighbor.get_voxel_pool()->species().serial() == loc)
+            {
+                tmp.push_back(neighbor);
+            }
+        }
+    }
+
     if (tmp.size() == 0)
     {
         return boost::none;
@@ -369,15 +380,9 @@ boost::optional<Voxel> SpatiocyteWorld::check_neighbor(const Voxel &voxel,
     return tmp[rng()->uniform_int(0, tmp.size() - 1)];
 }
 
-const Voxel SpatiocyteWorld::get_neighbor_randomly(const Voxel &voxel) const
-{
-    const Integer idx(rng()->uniform_int(0, num_neighbors(voxel) - 1));
-    return get_neighbor(voxel, idx);
-}
-
 const Voxel
 SpatiocyteWorld::get_neighbor_randomly(const Voxel &voxel,
-                                       Shape::dimension_kind dimension) const
+                                       Shape::dimension_kind dimension)
 {
     std::vector<Voxel> neighbors;
     for (Integer idx = 0; idx < num_neighbors(voxel); ++idx)
@@ -391,7 +396,15 @@ SpatiocyteWorld::get_neighbor_randomly(const Voxel &voxel,
     }
 
     const Integer idx(rng()->uniform_int(0, neighbors.size() - 1));
-    return neighbors.at(idx);
+    const auto neighbor = neighbors.at(idx);
+
+    if (const auto neighbors = interfaces_.find(neighbor))
+    {
+        const auto idx(rng()->uniform_int(0, neighbors->size() - 1));
+        return neighbors->at(idx);
+    }
+
+    return neighbor;
 }
 
 } // namespace spatiocyte

diff --git a/ecell4/spatiocyte/SpatiocyteWorld.hpp b/ecell4/spatiocyte/SpatiocyteWorld.hpp
@@ -767,9 +767,8 @@ class SpatiocyteWorld : public WorldInterface
                      voxel.space.lock()->get_neighbor(voxel.coordinate, nrand));
     }
 
-    const Voxel get_neighbor_randomly(const Voxel &voxel) const;
     const Voxel get_neighbor_randomly(const Voxel &voxel,
-                                      Shape::dimension_kind dimension) const;
+                                      Shape::dimension_kind dimension);
 
     const Species &draw_species(const Species &pttrn) const;
 

diff --git a/ecell4/spatiocyte/StepEvent.cpp b/ecell4/spatiocyte/StepEvent.cpp
@@ -67,7 +67,8 @@ void StepEvent3D::walk(const Real &alpha)
             continue;
         }
 
-        const Voxel neighbor(world_->get_neighbor_randomly(voxel));
+        const Voxel neighbor(
+            world_->get_neighbor_randomly(voxel, Shape::THREE));
 
         if (world_->can_move(voxel, neighbor))
         {
@@ -148,11 +149,6 @@ void StepEvent::attempt_reaction_(
     boost::shared_ptr<const VoxelPool> from_mt(voxel.get_voxel_pool());
     boost::shared_ptr<const VoxelPool> to_mt(dst.get_voxel_pool());
 
-    if (to_mt->is_vacant())
-    {
-        return;
-    }
-
     const Species &speciesA(from_mt->species());
     const Species &speciesB(to_mt->species());