From 177e20e23ed67492f675809ab336e9bdc44cab3f Mon Sep 17 00:00:00 2001
From: Suguru Kato <salmon.bioinformatics@gmail.com>
Date: Thu, 26 Dec 2019 16:06:52 +0900
Subject: [PATCH 01/21] :recycle: Add Voxel::get_neighbor_randomly()

---
 ecell4/spatiocyte/StepEvent.cpp | 3 +--
 ecell4/spatiocyte/Voxel.hpp     | 7 +++++++
 2 files changed, 8 insertions(+), 2 deletions(-)

diff --git a/ecell4/spatiocyte/StepEvent.cpp b/ecell4/spatiocyte/StepEvent.cpp
index 834020be4..7b07dc28d 100644
--- a/ecell4/spatiocyte/StepEvent.cpp
+++ b/ecell4/spatiocyte/StepEvent.cpp
@@ -49,7 +49,6 @@ void StepEvent3D::walk(const Real &alpha)
     for (const auto &info : voxels)
     {
         const Voxel voxel(world_->coordinate2voxel(info.coordinate));
-        const Integer rnd(rng->uniform_int(0, voxel.num_neighbors() - 1));
 
         if (voxel.get_voxel_pool() != mpool_)
         {
@@ -58,7 +57,7 @@ void StepEvent3D::walk(const Real &alpha)
             continue;
         }
 
-        const Voxel neighbor(voxel.get_neighbor(rnd));
+        const Voxel neighbor(voxel.get_neighbor_randomly(rng));
 
         if (world_->can_move(voxel, neighbor))
         {
diff --git a/ecell4/spatiocyte/Voxel.hpp b/ecell4/spatiocyte/Voxel.hpp
index b75f1ffd6..3081f482d 100644
--- a/ecell4/spatiocyte/Voxel.hpp
+++ b/ecell4/spatiocyte/Voxel.hpp
@@ -45,6 +45,13 @@ struct Voxel
     {
         return Voxel(space, space.lock()->get_neighbor(coordinate, nrand));
     }
+
+    Voxel get_neighbor_randomly(
+        const boost::shared_ptr<RandomNumberGenerator> &rng) const
+    {
+        const Integer idx(rng->uniform_int(0, num_neighbors() - 1));
+        return Voxel(space, space.lock()->get_neighbor(coordinate, idx));
+    }
 };
 
 } // namespace spatiocyte

From 75772086ae7b34732e4496ba85e8b3c59d03112c Mon Sep 17 00:00:00 2001
From: Suguru Kato <salmon.bioinformatics@gmail.com>
Date: Thu, 26 Dec 2019 17:12:44 +0900
Subject: [PATCH 02/21] :recycle: Add another Voxel::get_neighbor_randomly()

For StepEvent2D
---
 ecell4/spatiocyte/SpatiocyteWorld.hpp | 26 +++++++++++++++++-----
 ecell4/spatiocyte/Voxel.cpp           | 32 +++++++++++++++++++++++++++
 ecell4/spatiocyte/Voxel.hpp           | 17 ++++++++++----
 3 files changed, 66 insertions(+), 9 deletions(-)
 create mode 100644 ecell4/spatiocyte/Voxel.cpp

diff --git a/ecell4/spatiocyte/SpatiocyteWorld.hpp b/ecell4/spatiocyte/SpatiocyteWorld.hpp
index 10b42626c..bfd18deb9 100644
--- a/ecell4/spatiocyte/SpatiocyteWorld.hpp
+++ b/ecell4/spatiocyte/SpatiocyteWorld.hpp
@@ -35,6 +35,7 @@ struct MoleculeInfo
     const std::string loc;
     const Shape::dimension_kind dimension;
 };
+
 class SpatiocyteWorld : public WorldInterface
 {
 public:
@@ -240,7 +241,7 @@ class SpatiocyteWorld : public WorldInterface
         for (const auto &space : spaces_)
         {
             if (const auto coordinate = space->get_coordinate(pid))
-                return Voxel(space, *coordinate);
+                return Voxel(this, space, *coordinate);
         }
         return boost::none;
     }
@@ -443,7 +444,7 @@ class SpatiocyteWorld : public WorldInterface
      */
     Voxel get_voxel_nearby(const Real3 &pos) const
     {
-        return Voxel(get_root(), get_root()->position2coordinate(pos));
+        return Voxel(this, get_root(), get_root()->position2coordinate(pos));
     }
 
     /*
@@ -513,6 +514,16 @@ class SpatiocyteWorld : public WorldInterface
      * SpatiocyteWorld API
      */
 
+    const MoleculeInfo get_molecule_info(const Species &sp) const
+    {
+        const auto itr = molecule_info_cache_.find(sp);
+        if (itr != molecule_info_cache_.end())
+        {
+            return itr->second;
+        }
+        throw NotFound("MoleculeInfo not found");
+    }
+
     /**
      * draw attributes of species and return it as a molecule info.
      * @param sp a species
@@ -775,6 +786,11 @@ class SpatiocyteWorld : public WorldInterface
         return get_molecule_info(species).dimension;
     }
 
+    Shape::dimension_kind get_dimension(const Species &species) const
+    {
+        return get_molecule_info(species).dimension;
+    }
+
     /**
      * static members
      */
@@ -879,9 +895,9 @@ class SpatiocyteWorld : public WorldInterface
     std::vector<ParticleBase<Voxel>> list_voxels_private(ListFn list_fn) const
     {
         return map_voxels<ParticleBase<Voxel>>(
-            list_fn, [](const space_type &space, const VoxelView &view) {
+            list_fn, [this](const space_type &space, const VoxelView &view) {
                 return ParticleBase<Voxel>(view.pid, view.species,
-                                           Voxel(space, view.voxel));
+                                           Voxel(this, space, view.voxel));
             });
     }
 
@@ -890,7 +906,7 @@ class SpatiocyteWorld : public WorldInterface
     Voxel coordinate2voxel(const coordinate_type &coordinate)
     {
         coordinate_type coord(coordinate);
-        return Voxel(get_space(coord), coord);
+        return Voxel(this, get_space(coord), coord);
     }
 
 protected:
diff --git a/ecell4/spatiocyte/Voxel.cpp b/ecell4/spatiocyte/Voxel.cpp
new file mode 100644
index 000000000..dd3237927
--- /dev/null
+++ b/ecell4/spatiocyte/Voxel.cpp
@@ -0,0 +1,32 @@
+#include "SpatiocyteWorld.hpp"
+#include "Voxel.hpp"
+
+namespace ecell4
+{
+
+namespace spatiocyte
+{
+
+Voxel Voxel::get_neighbor_randomly(
+    const boost::shared_ptr<RandomNumberGenerator> &rng,
+    Shape::dimension_kind dimension) const
+{
+    std::vector<Voxel> neighbors;
+    for (Integer idx = 0; idx < num_neighbors(); ++idx)
+    {
+        const Voxel neighbor = get_neighbor(idx);
+        if (world->get_dimension(neighbor.get_voxel_pool()->species()) >
+            dimension)
+        {
+            continue;
+        }
+        neighbors.push_back(neighbor);
+    }
+
+    const Integer idx(rng->uniform_int(0, neighbors.size() - 1));
+    return neighbors.at(idx);
+}
+
+} // namespace spatiocyte
+
+} // namespace ecell4
diff --git a/ecell4/spatiocyte/Voxel.hpp b/ecell4/spatiocyte/Voxel.hpp
index 3081f482d..861cd580b 100644
--- a/ecell4/spatiocyte/Voxel.hpp
+++ b/ecell4/spatiocyte/Voxel.hpp
@@ -11,15 +11,19 @@ namespace ecell4
 namespace spatiocyte
 {
 
+class SpatiocyteWorld;
+
 struct Voxel
 {
     typedef Integer coordinate_type;
 
-    Voxel(boost::weak_ptr<VoxelSpaceBase> space, coordinate_type coordinate)
-        : space(space), coordinate(coordinate)
+    Voxel(const SpatiocyteWorld *world, boost::weak_ptr<VoxelSpaceBase> space,
+          coordinate_type coordinate)
+        : world(world), space(space), coordinate(coordinate)
     {
     }
 
+    const SpatiocyteWorld *world;
     boost::weak_ptr<VoxelSpaceBase> space;
     coordinate_type coordinate;
 
@@ -43,15 +47,20 @@ struct Voxel
 
     Voxel get_neighbor(Integer nrand) const
     {
-        return Voxel(space, space.lock()->get_neighbor(coordinate, nrand));
+        return Voxel(world, space,
+                     space.lock()->get_neighbor(coordinate, nrand));
     }
 
     Voxel get_neighbor_randomly(
         const boost::shared_ptr<RandomNumberGenerator> &rng) const
     {
         const Integer idx(rng->uniform_int(0, num_neighbors() - 1));
-        return Voxel(space, space.lock()->get_neighbor(coordinate, idx));
+        return get_neighbor(idx);
     }
+
+    Voxel
+    get_neighbor_randomly(const boost::shared_ptr<RandomNumberGenerator> &rng,
+                          Shape::dimension_kind dimension) const;
 };
 
 } // namespace spatiocyte

From 2dc40f8716b6a68cbde94a9afff2379f4deb3610 Mon Sep 17 00:00:00 2001
From: Suguru Kato <salmon.bioinformatics@gmail.com>
Date: Thu, 26 Dec 2019 18:46:27 +0900
Subject: [PATCH 03/21] :recycle: StepEvent2D uses
 Voxel::get_neighbor_randomly()

---
 ecell4/spatiocyte/SpatiocyteEvent.hpp |  3 ---
 ecell4/spatiocyte/StepEvent.cpp       | 35 +++++++--------------------
 2 files changed, 9 insertions(+), 29 deletions(-)

diff --git a/ecell4/spatiocyte/SpatiocyteEvent.hpp b/ecell4/spatiocyte/SpatiocyteEvent.hpp
index c51b402f4..faf498c7a 100644
--- a/ecell4/spatiocyte/SpatiocyteEvent.hpp
+++ b/ecell4/spatiocyte/SpatiocyteEvent.hpp
@@ -87,9 +87,6 @@ struct StepEvent2D : StepEvent
                 const Species &species, const Real &t, const Real alpha = 1.0);
 
     void walk(const Real &alpha);
-
-protected:
-    std::vector<unsigned int> nids_; // neighbor indexes
 };
 
 struct ZerothOrderReactionEvent : SpatiocyteEvent
diff --git a/ecell4/spatiocyte/StepEvent.cpp b/ecell4/spatiocyte/StepEvent.cpp
index 7b07dc28d..281693506 100644
--- a/ecell4/spatiocyte/StepEvent.cpp
+++ b/ecell4/spatiocyte/StepEvent.cpp
@@ -89,10 +89,6 @@ StepEvent2D::StepEvent2D(boost::shared_ptr<Model> model,
         dt_ = R * R / D * alpha_;
 
     time_ = t + dt_;
-
-    nids_.clear();
-    for (unsigned int i(0); i < 12; ++i)
-        nids_.push_back(i);
 }
 
 void StepEvent2D::walk(const Real &alpha)
@@ -120,31 +116,18 @@ void StepEvent2D::walk(const Real &alpha)
         }
 
         const std::size_t num_neighbors(voxel.num_neighbors());
+        const Voxel neighbor(voxel.get_neighbor_randomly(rng, Shape::TWO));
 
-        ecell4::shuffle(*(rng.get()), nids_);
-        for (const auto &index : nids_)
+        if (world_->can_move(voxel, neighbor))
         {
-            if (index >= num_neighbors)
-                continue;
-
-            const Voxel neighbor(voxel.get_neighbor(index));
-            boost::shared_ptr<const VoxelPool> target(
-                neighbor.get_voxel_pool());
-
-            if (world_->get_dimension(target->species()) > Shape::TWO)
-                continue;
-
-            if (world_->can_move(voxel, neighbor))
-            {
-                if (rng->uniform(0, 1) <= alpha)
-                    world_->move(voxel, neighbor, /*candidate=*/idx);
-            }
-            else
-            {
-                attempt_reaction_(info, neighbor, alpha);
-            }
-            break;
+            if (rng->uniform(0, 1) <= alpha)
+                world_->move(voxel, neighbor, /*candidate=*/idx);
+        }
+        else
+        {
+            attempt_reaction_(info, neighbor, alpha);
         }
+
         ++idx;
     }
 }

From 2eb17252c69be521f13d6618aa54eb2ee946824b Mon Sep 17 00:00:00 2001
From: Suguru Kato <salmon.bioinformatics@gmail.com>
Date: Thu, 9 Jan 2020 10:31:45 +0900
Subject: [PATCH 04/21] :recycle: Remove rng from args of
 get_neighbor_randomly()

---
 ecell4/spatiocyte/SpatiocyteWorld.hpp |  2 +-
 ecell4/spatiocyte/StepEvent.cpp       |  4 ++--
 ecell4/spatiocyte/Voxel.cpp           | 12 ++++++++----
 ecell4/spatiocyte/Voxel.hpp           | 11 ++---------
 4 files changed, 13 insertions(+), 16 deletions(-)

diff --git a/ecell4/spatiocyte/SpatiocyteWorld.hpp b/ecell4/spatiocyte/SpatiocyteWorld.hpp
index bfd18deb9..9dc180d1f 100644
--- a/ecell4/spatiocyte/SpatiocyteWorld.hpp
+++ b/ecell4/spatiocyte/SpatiocyteWorld.hpp
@@ -750,7 +750,7 @@ class SpatiocyteWorld : public WorldInterface
 
     const Species &draw_species(const Species &pttrn) const;
 
-    boost::shared_ptr<RandomNumberGenerator> rng() { return rng_; }
+    boost::shared_ptr<RandomNumberGenerator> rng() const { return rng_; }
 
     void bind_to(boost::shared_ptr<Model> model)
     {
diff --git a/ecell4/spatiocyte/StepEvent.cpp b/ecell4/spatiocyte/StepEvent.cpp
index 281693506..44df254d6 100644
--- a/ecell4/spatiocyte/StepEvent.cpp
+++ b/ecell4/spatiocyte/StepEvent.cpp
@@ -57,7 +57,7 @@ void StepEvent3D::walk(const Real &alpha)
             continue;
         }
 
-        const Voxel neighbor(voxel.get_neighbor_randomly(rng));
+        const Voxel neighbor(voxel.get_neighbor_randomly());
 
         if (world_->can_move(voxel, neighbor))
         {
@@ -116,7 +116,7 @@ void StepEvent2D::walk(const Real &alpha)
         }
 
         const std::size_t num_neighbors(voxel.num_neighbors());
-        const Voxel neighbor(voxel.get_neighbor_randomly(rng, Shape::TWO));
+        const Voxel neighbor(voxel.get_neighbor_randomly(Shape::TWO));
 
         if (world_->can_move(voxel, neighbor))
         {
diff --git a/ecell4/spatiocyte/Voxel.cpp b/ecell4/spatiocyte/Voxel.cpp
index dd3237927..3c091c911 100644
--- a/ecell4/spatiocyte/Voxel.cpp
+++ b/ecell4/spatiocyte/Voxel.cpp
@@ -7,9 +7,13 @@ namespace ecell4
 namespace spatiocyte
 {
 
-Voxel Voxel::get_neighbor_randomly(
-    const boost::shared_ptr<RandomNumberGenerator> &rng,
-    Shape::dimension_kind dimension) const
+Voxel Voxel::get_neighbor_randomly() const
+{
+    const Integer idx(world->rng()->uniform_int(0, num_neighbors() - 1));
+    return get_neighbor(idx);
+}
+
+Voxel Voxel::get_neighbor_randomly(Shape::dimension_kind dimension) const
 {
     std::vector<Voxel> neighbors;
     for (Integer idx = 0; idx < num_neighbors(); ++idx)
@@ -23,7 +27,7 @@ Voxel Voxel::get_neighbor_randomly(
         neighbors.push_back(neighbor);
     }
 
-    const Integer idx(rng->uniform_int(0, neighbors.size() - 1));
+    const Integer idx(world->rng()->uniform_int(0, neighbors.size() - 1));
     return neighbors.at(idx);
 }
 
diff --git a/ecell4/spatiocyte/Voxel.hpp b/ecell4/spatiocyte/Voxel.hpp
index 861cd580b..9e3406ad0 100644
--- a/ecell4/spatiocyte/Voxel.hpp
+++ b/ecell4/spatiocyte/Voxel.hpp
@@ -51,16 +51,9 @@ struct Voxel
                      space.lock()->get_neighbor(coordinate, nrand));
     }
 
-    Voxel get_neighbor_randomly(
-        const boost::shared_ptr<RandomNumberGenerator> &rng) const
-    {
-        const Integer idx(rng->uniform_int(0, num_neighbors() - 1));
-        return get_neighbor(idx);
-    }
+    Voxel get_neighbor_randomly() const;
 
-    Voxel
-    get_neighbor_randomly(const boost::shared_ptr<RandomNumberGenerator> &rng,
-                          Shape::dimension_kind dimension) const;
+    Voxel get_neighbor_randomly(Shape::dimension_kind dimension) const;
 };
 
 } // namespace spatiocyte

From 41807bbbde0649354f177bb28cc2fdcec215205b Mon Sep 17 00:00:00 2001
From: Suguru Kato <salmon.bioinformatics@gmail.com>
Date: Thu, 16 Jan 2020 11:13:37 +0900
Subject: [PATCH 05/21] :recycle: Move functions for neighbor to
 SpatiocyteWorld

---
 ecell4/python_api/spatiocyte.cpp              | 175 +++++++++---------
 ecell4/spatiocyte/SpatiocyteWorld.cpp         |  37 +++-
 ecell4/spatiocyte/SpatiocyteWorld.hpp         |  23 ++-
 ecell4/spatiocyte/StepEvent.cpp               |   5 +-
 ecell4/spatiocyte/Voxel.cpp                   |  36 ----
 ecell4/spatiocyte/Voxel.hpp                   |  21 +--
 .../spatiocyte/tests/SpatiocyteWorld_test.cpp |   4 +-
 7 files changed, 147 insertions(+), 154 deletions(-)
 delete mode 100644 ecell4/spatiocyte/Voxel.cpp

diff --git a/ecell4/python_api/spatiocyte.cpp b/ecell4/python_api/spatiocyte.cpp
index df0ff028f..c73318584 100644
--- a/ecell4/python_api/spatiocyte.cpp
+++ b/ecell4/python_api/spatiocyte.cpp
@@ -1,12 +1,12 @@
 #include "python_api.hpp"
 
 #include <ecell4/core/OffLatticeSpace.hpp>
-#include <ecell4/spatiocyte/SpatiocyteReactions.hpp>
+#include <ecell4/spatiocyte/OffLattice.hpp>
 #include <ecell4/spatiocyte/SpatiocyteFactory.hpp>
+#include <ecell4/spatiocyte/SpatiocyteReactions.hpp>
 #include <ecell4/spatiocyte/SpatiocyteSimulator.hpp>
 #include <ecell4/spatiocyte/SpatiocyteWorld.hpp>
 #include <ecell4/spatiocyte/Voxel.hpp>
-#include <ecell4/spatiocyte/OffLattice.hpp>
 
 #include "simulator.hpp"
 #include "simulator_factory.hpp"
@@ -21,32 +21,28 @@ namespace ecell4
 namespace python_api
 {
 
-static inline
-void define_reaction_info(py::module& m)
+static inline void define_reaction_info(py::module &m)
 {
     using container_type = ReactionInfo::container_type;
     py::class_<ReactionInfo>(m, "ReactionInfo")
-        .def(py::init<const Real, const container_type&, const container_type&>(),
-                py::arg("t"), py::arg("reactants"), py::arg("products"))
+        .def(py::init<const Real, const container_type &,
+                      const container_type &>(),
+             py::arg("t"), py::arg("reactants"), py::arg("products"))
         .def("t", &ReactionInfo::t)
         .def("reactants", &ReactionInfo::reactants)
         .def("products", &ReactionInfo::products)
         .def(py::pickle(
-            [](const ReactionInfo& self)
-            {
-                return py::make_tuple(self.t(), self.reactants(), self.products());
+            [](const ReactionInfo &self) {
+                return py::make_tuple(self.t(), self.reactants(),
+                                      self.products());
             },
-            [](py::tuple t)
-            {
+            [](py::tuple t) {
                 if (t.size() != 3)
                     throw std::runtime_error("Invalid state");
-                return ReactionInfo(
-                    t[0].cast<Real>(),
-                    t[1].cast<container_type>(),
-                    t[2].cast<container_type>()
-                );
-            }
-        ));
+                return ReactionInfo(t[0].cast<Real>(),
+                                    t[1].cast<container_type>(),
+                                    t[2].cast<container_type>());
+            }));
 
     py::class_<ReactionInfo::Item>(m, "ReactionInfoItem")
         .def_readonly("pid", &ReactionInfo::Item::pid)
@@ -54,28 +50,28 @@ void define_reaction_info(py::module& m)
         .def_readonly("voxel", &ReactionInfo::Item::voxel);
 }
 
-static inline
-void define_spatiocyte_factory(py::module& m)
+static inline void define_spatiocyte_factory(py::module &m)
 {
     py::class_<SpatiocyteFactory> factory(m, "SpatiocyteFactory");
     factory
         .def(py::init<const Real>(),
-                py::arg("voxel_radius") = SpatiocyteFactory::default_voxel_radius())
+             py::arg("voxel_radius") =
+                 SpatiocyteFactory::default_voxel_radius())
         .def("rng", &SpatiocyteFactory::rng);
     define_factory_functions(factory);
 
     m.attr("Factory") = factory;
 }
 
-static inline
-void define_spatiocyte_simulator(py::module& m)
+static inline void define_spatiocyte_simulator(py::module &m)
 {
     py::class_<SpatiocyteSimulator, Simulator, PySimulator<SpatiocyteSimulator>,
-        boost::shared_ptr<SpatiocyteSimulator>> simulator(m, "SpatiocyteSimulator");
-    simulator
-        .def(py::init<boost::shared_ptr<SpatiocyteWorld>>(), py::arg("w"))
-        .def(py::init<boost::shared_ptr<SpatiocyteWorld>, boost::shared_ptr<Model>>(),
-                py::arg("w"), py::arg("m"))
+               boost::shared_ptr<SpatiocyteSimulator>>
+        simulator(m, "SpatiocyteSimulator");
+    simulator.def(py::init<boost::shared_ptr<SpatiocyteWorld>>(), py::arg("w"))
+        .def(py::init<boost::shared_ptr<SpatiocyteWorld>,
+                      boost::shared_ptr<Model>>(),
+             py::arg("w"), py::arg("m"))
         .def("last_reactions", &SpatiocyteSimulator::last_reactions)
         .def("set_t", &SpatiocyteSimulator::set_t);
     define_simulator_functions(simulator);
@@ -83,34 +79,40 @@ void define_spatiocyte_simulator(py::module& m)
     m.attr("Simulator") = simulator;
 }
 
-static inline
-void define_spatiocyte_world(py::module& m)
+static inline void define_spatiocyte_world(py::module &m)
 {
-    py::class_<SpatiocyteWorld, ecell4::WorldInterface, PyWorldImpl<SpatiocyteWorld>,
-        boost::shared_ptr<SpatiocyteWorld>> world(m, "SpatiocyteWorld");
+    py::class_<SpatiocyteWorld, ecell4::WorldInterface,
+               PyWorldImpl<SpatiocyteWorld>, boost::shared_ptr<SpatiocyteWorld>>
+        world(m, "SpatiocyteWorld");
     world
-        .def(py::init<const Real3&>(),
-                py::arg("edge_lengths") = Real3(1.0, 1.0, 1.0))
-        .def(py::init<const Real3&, const Real&>(),
-                py::arg("edge_lengths"), py::arg("voxel_radius"))
-        .def(py::init<const Real3&, const Real&, const boost::shared_ptr<RandomNumberGenerator>&>(),
-                py::arg("edge_lengths"),
-                py::arg("voxel_radius"),
-                py::arg("rng"))
+        .def(py::init<const Real3 &>(),
+             py::arg("edge_lengths") = Real3(1.0, 1.0, 1.0))
+        .def(py::init<const Real3 &, const Real &>(), py::arg("edge_lengths"),
+             py::arg("voxel_radius"))
+        .def(py::init<const Real3 &, const Real &,
+                      const boost::shared_ptr<RandomNumberGenerator> &>(),
+             py::arg("edge_lengths"), py::arg("voxel_radius"), py::arg("rng"))
         .def(py::init<const std::string>(), py::arg("filename"))
 
-        .def("new_particle", (boost::optional<ParticleID> (SpatiocyteWorld::*)(const Particle&)) &SpatiocyteWorld::new_particle)
-        .def("new_particle", (boost::optional<ParticleID> (SpatiocyteWorld::*)(const Species&, const Real3&)) &SpatiocyteWorld::new_particle)
+        .def("new_particle", (boost::optional<ParticleID>(SpatiocyteWorld::*)(
+                                 const Particle &)) &
+                                 SpatiocyteWorld::new_particle)
+        .def("new_particle", (boost::optional<ParticleID>(SpatiocyteWorld::*)(
+                                 const Species &, const Real3 &)) &
+                                 SpatiocyteWorld::new_particle)
         .def("remove_particle", &SpatiocyteWorld::remove_particle)
-        .def("list_structure_particles", &SpatiocyteWorld::list_structure_particles)
-        .def("list_non_structure_particles", &SpatiocyteWorld::list_non_structure_particles)
+        .def("list_structure_particles",
+             &SpatiocyteWorld::list_structure_particles)
+        .def("list_non_structure_particles",
+             &SpatiocyteWorld::list_non_structure_particles)
         .def("update_particle", &SpatiocyteWorld::update_particle)
         .def("add_molecules",
-                (bool (SpatiocyteWorld::*)(const Species&, const Integer&))
-                &SpatiocyteWorld::add_molecules)
+             (bool (SpatiocyteWorld::*)(const Species &, const Integer &)) &
+                 SpatiocyteWorld::add_molecules)
         .def("add_molecules",
-                (bool (SpatiocyteWorld::*)(const Species&, const Integer&, const boost::shared_ptr<const Shape>))
-                &SpatiocyteWorld::add_molecules)
+             (bool (SpatiocyteWorld::*)(const Species &, const Integer &,
+                                        const boost::shared_ptr<const Shape>)) &
+                 SpatiocyteWorld::add_molecules)
         .def("remove_molecules", &SpatiocyteWorld::remove_molecules)
         .def("voxel_volume", &SpatiocyteWorld::voxel_volume)
         .def("get_volume", &SpatiocyteWorld::get_volume)
@@ -118,14 +120,17 @@ void define_spatiocyte_world(py::module& m)
         .def("get_species_at", &SpatiocyteWorld::get_species_at)
         .def("has_particle_at", &SpatiocyteWorld::has_particle_at)
         .def("set_value", &SpatiocyteWorld::set_value)
-        .def("new_particle", (boost::optional<ParticleID> (SpatiocyteWorld::*)(const Species&, const Voxel&))&SpatiocyteWorld::new_particle)
+        .def("new_particle", (boost::optional<ParticleID>(SpatiocyteWorld::*)(
+                                 const Species &, const Voxel &)) &
+                                 SpatiocyteWorld::new_particle)
         .def("new_voxel",
-            (boost::optional<ParticleID> (SpatiocyteWorld::*)(const Species&, const Voxel&))&SpatiocyteWorld::new_particle,
-            R"pbdoc(
+             (boost::optional<ParticleID>(SpatiocyteWorld::*)(const Species &,
+                                                              const Voxel &)) &
+                 SpatiocyteWorld::new_particle,
+             R"pbdoc(
                 .. deprecated::3.0
                    Use :func:`new_particle` instead.
-            )pbdoc"
-        )
+            )pbdoc")
         .def("new_voxel_structure", &SpatiocyteWorld::new_voxel_structure)
         .def("voxel_radius", &SpatiocyteWorld::voxel_radius)
         .def("size", &SpatiocyteWorld::size)
@@ -148,54 +153,57 @@ void define_spatiocyte_world(py::module& m)
         )pbdoc")
         .def("rng", &SpatiocyteWorld::rng)
         .def("add_offlattice",
-            [](SpatiocyteWorld& self, const Species& species, const OffLattice& offlattice)
-            {
-                self.add_space(offlattice.generate_space(species));
-            })
-        .def_static("calculate_voxel_volume", &SpatiocyteWorld::calculate_voxel_volume)
-        .def_static("calculate_hcp_lengths", &SpatiocyteWorld::calculate_hcp_lengths)
+             [](SpatiocyteWorld &self, const Species &species,
+                const OffLattice &offlattice) {
+                 self.add_space(offlattice.generate_space(species));
+             })
+        .def_static("calculate_voxel_volume",
+                    &SpatiocyteWorld::calculate_voxel_volume)
+        .def_static("calculate_hcp_lengths",
+                    &SpatiocyteWorld::calculate_hcp_lengths)
         .def_static("calculate_shape", &SpatiocyteWorld::calculate_shape)
-        .def_static("calculate_volume", &SpatiocyteWorld::calculate_volume);
-
-    m.def("create_spatiocyte_world_cell_list_impl", &create_spatiocyte_world_cell_list_impl);
-    m.def("create_spatiocyte_world_vector_impl", &create_spatiocyte_world_vector_impl);
-    m.def("create_spatiocyte_world_square_offlattice_impl", &allocate_spatiocyte_world_square_offlattice_impl);
+        .def_static("calculate_volume", &SpatiocyteWorld::calculate_volume)
+        .def("list_neighbors",
+             [](const SpatiocyteWorld &self, const Voxel &voxel) {
+                 std::vector<Voxel> list;
+                 for (auto i = 0; i < self.num_neighbors(voxel); ++i)
+                 {
+                     list.push_back(self.get_neighbor(voxel, i));
+                 }
+                 return list;
+             });
+
+    m.def("create_spatiocyte_world_cell_list_impl",
+          &create_spatiocyte_world_cell_list_impl);
+    m.def("create_spatiocyte_world_vector_impl",
+          &create_spatiocyte_world_vector_impl);
+    m.def("create_spatiocyte_world_square_offlattice_impl",
+          &allocate_spatiocyte_world_square_offlattice_impl);
 
     m.attr("World") = world;
 }
 
-static inline
-void define_voxel(py::module& m)
+static inline void define_voxel(py::module &m)
 {
     py::class_<Voxel>(m, "Voxel")
-        .def("position", &Voxel::position)
 #ifndef NDEBUG
         .def_readonly("coordinate", &Voxel::coordinate)
 #endif
-        .def("list_neighbors",
-            [](const Voxel& self)
-            {
-                std::vector<Voxel> list;
-                for (auto i = 0; i < self.num_neighbors(); ++i)
-                {
-                    list.push_back(self.get_neighbor(i));
-                }
-                return list;
-            });
+        .def("position", &Voxel::position);
 }
 
-static inline
-void define_offlattice(py::module& m)
+static inline void define_offlattice(py::module &m)
 {
     py::class_<OffLattice>(m, "OffLattice")
-        .def(py::init<const Real, const OffLattice::positions_type, const OffLattice::adjoining_pairs_type>())
+        .def(py::init<const Real, const OffLattice::positions_type,
+                      const OffLattice::adjoining_pairs_type>())
         .def(py::init<const Real, const OffLattice::positions_type>())
         .def("voxel_radius", &OffLattice::voxel_radius)
         .def("positions", &OffLattice::positions)
         .def("adjoining_pairs", &OffLattice::adjoining_pairs);
 }
 
-void setup_spatiocyte_module(py::module& m)
+void setup_spatiocyte_module(py::module &m)
 {
     define_reaction_info(m);
     define_offlattice(m);
@@ -205,7 +213,6 @@ void setup_spatiocyte_module(py::module& m)
     define_voxel(m);
 }
 
-}
-
-}
+} // namespace python_api
 
+} // namespace ecell4
diff --git a/ecell4/spatiocyte/SpatiocyteWorld.cpp b/ecell4/spatiocyte/SpatiocyteWorld.cpp
index a29b7305a..56a423d95 100644
--- a/ecell4/spatiocyte/SpatiocyteWorld.cpp
+++ b/ecell4/spatiocyte/SpatiocyteWorld.cpp
@@ -286,8 +286,8 @@ bool SpatiocyteWorld::is_surface_voxel(
     if (L > 0 || L < -2 * voxel_radius())
         return false;
 
-    for (Integer i(0); i < voxel.num_neighbors(); ++i)
-        if (shape->is_inside(voxel.get_neighbor(i).position()) > 0)
+    for (Integer i(0); i < num_neighbors(voxel); ++i)
+        if (shape->is_inside(get_neighbor(voxel, i).position()) > 0)
             return true;
 
     return false;
@@ -322,14 +322,14 @@ void SpatiocyteWorld::remove_molecules(const Species &sp, const Integer &num)
 boost::optional<Voxel> SpatiocyteWorld::check_neighbor(const Voxel &voxel,
                                                        const std::string &loc)
 {
-    const std::size_t num_neighbors(voxel.num_neighbors());
+    const std::size_t num(num_neighbors(voxel));
 
     std::vector<Voxel> tmp;
-    tmp.reserve(num_neighbors);
+    tmp.reserve(num);
 
-    for (unsigned int rnd(0); rnd < num_neighbors; ++rnd)
+    for (unsigned int rnd(0); rnd < num; ++rnd)
     {
-        const Voxel neighbor(voxel.get_neighbor(rnd));
+        const Voxel neighbor(get_neighbor(voxel, rnd));
         boost::shared_ptr<const VoxelPool> mt(neighbor.get_voxel_pool());
         const std::string serial(mt->is_vacant() ? "" : mt->species().serial());
         if (serial == loc)
@@ -346,6 +346,31 @@ 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
+{
+    std::vector<Voxel> neighbors;
+    for (Integer idx = 0; idx < num_neighbors(voxel); ++idx)
+    {
+        const Voxel neighbor = get_neighbor(voxel, idx);
+        if (get_dimension(neighbor.get_voxel_pool()->species()) > dimension)
+        {
+            continue;
+        }
+        neighbors.push_back(neighbor);
+    }
+
+    const Integer idx(rng()->uniform_int(0, neighbors.size() - 1));
+    return neighbors.at(idx);
+}
+
 } // namespace spatiocyte
 
 } // namespace ecell4
diff --git a/ecell4/spatiocyte/SpatiocyteWorld.hpp b/ecell4/spatiocyte/SpatiocyteWorld.hpp
index 9dc180d1f..25eac5432 100644
--- a/ecell4/spatiocyte/SpatiocyteWorld.hpp
+++ b/ecell4/spatiocyte/SpatiocyteWorld.hpp
@@ -241,7 +241,7 @@ class SpatiocyteWorld : public WorldInterface
         for (const auto &space : spaces_)
         {
             if (const auto coordinate = space->get_coordinate(pid))
-                return Voxel(this, space, *coordinate);
+                return Voxel(space, *coordinate);
         }
         return boost::none;
     }
@@ -444,7 +444,7 @@ class SpatiocyteWorld : public WorldInterface
      */
     Voxel get_voxel_nearby(const Real3 &pos) const
     {
-        return Voxel(this, get_root(), get_root()->position2coordinate(pos));
+        return Voxel(get_root(), get_root()->position2coordinate(pos));
     }
 
     /*
@@ -746,6 +746,21 @@ class SpatiocyteWorld : public WorldInterface
     boost::optional<Voxel> check_neighbor(const Voxel &voxel,
                                           const std::string &loc);
 
+    const Integer num_neighbors(const Voxel &voxel) const
+    {
+        return voxel.space.lock()->num_neighbors(voxel.coordinate);
+    }
+
+    const Voxel get_neighbor(const Voxel &voxel, Integer nrand) const
+    {
+        return Voxel(voxel.space,
+                     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;
+
     // bool update_molecule(coordinate_type at, Species species);
 
     const Species &draw_species(const Species &pttrn) const;
@@ -897,7 +912,7 @@ class SpatiocyteWorld : public WorldInterface
         return map_voxels<ParticleBase<Voxel>>(
             list_fn, [this](const space_type &space, const VoxelView &view) {
                 return ParticleBase<Voxel>(view.pid, view.species,
-                                           Voxel(this, space, view.voxel));
+                                           Voxel(space, view.voxel));
             });
     }
 
@@ -906,7 +921,7 @@ class SpatiocyteWorld : public WorldInterface
     Voxel coordinate2voxel(const coordinate_type &coordinate)
     {
         coordinate_type coord(coordinate);
-        return Voxel(this, get_space(coord), coord);
+        return Voxel(get_space(coord), coord);
     }
 
 protected:
diff --git a/ecell4/spatiocyte/StepEvent.cpp b/ecell4/spatiocyte/StepEvent.cpp
index 44df254d6..f5cad1e04 100644
--- a/ecell4/spatiocyte/StepEvent.cpp
+++ b/ecell4/spatiocyte/StepEvent.cpp
@@ -57,7 +57,7 @@ void StepEvent3D::walk(const Real &alpha)
             continue;
         }
 
-        const Voxel neighbor(voxel.get_neighbor_randomly());
+        const Voxel neighbor(world_->get_neighbor_randomly(voxel));
 
         if (world_->can_move(voxel, neighbor))
         {
@@ -115,8 +115,7 @@ void StepEvent2D::walk(const Real &alpha)
             continue;
         }
 
-        const std::size_t num_neighbors(voxel.num_neighbors());
-        const Voxel neighbor(voxel.get_neighbor_randomly(Shape::TWO));
+        const Voxel neighbor(world_->get_neighbor_randomly(voxel, Shape::TWO));
 
         if (world_->can_move(voxel, neighbor))
         {
diff --git a/ecell4/spatiocyte/Voxel.cpp b/ecell4/spatiocyte/Voxel.cpp
deleted file mode 100644
index 3c091c911..000000000
--- a/ecell4/spatiocyte/Voxel.cpp
+++ /dev/null
@@ -1,36 +0,0 @@
-#include "SpatiocyteWorld.hpp"
-#include "Voxel.hpp"
-
-namespace ecell4
-{
-
-namespace spatiocyte
-{
-
-Voxel Voxel::get_neighbor_randomly() const
-{
-    const Integer idx(world->rng()->uniform_int(0, num_neighbors() - 1));
-    return get_neighbor(idx);
-}
-
-Voxel Voxel::get_neighbor_randomly(Shape::dimension_kind dimension) const
-{
-    std::vector<Voxel> neighbors;
-    for (Integer idx = 0; idx < num_neighbors(); ++idx)
-    {
-        const Voxel neighbor = get_neighbor(idx);
-        if (world->get_dimension(neighbor.get_voxel_pool()->species()) >
-            dimension)
-        {
-            continue;
-        }
-        neighbors.push_back(neighbor);
-    }
-
-    const Integer idx(world->rng()->uniform_int(0, neighbors.size() - 1));
-    return neighbors.at(idx);
-}
-
-} // namespace spatiocyte
-
-} // namespace ecell4
diff --git a/ecell4/spatiocyte/Voxel.hpp b/ecell4/spatiocyte/Voxel.hpp
index 9e3406ad0..938997272 100644
--- a/ecell4/spatiocyte/Voxel.hpp
+++ b/ecell4/spatiocyte/Voxel.hpp
@@ -17,13 +17,11 @@ struct Voxel
 {
     typedef Integer coordinate_type;
 
-    Voxel(const SpatiocyteWorld *world, boost::weak_ptr<VoxelSpaceBase> space,
-          coordinate_type coordinate)
-        : world(world), space(space), coordinate(coordinate)
+    Voxel(boost::weak_ptr<VoxelSpaceBase> space, coordinate_type coordinate)
+        : space(space), coordinate(coordinate)
     {
     }
 
-    const SpatiocyteWorld *world;
     boost::weak_ptr<VoxelSpaceBase> space;
     coordinate_type coordinate;
 
@@ -39,21 +37,6 @@ struct Voxel
     {
         return space.lock()->get_voxel_pool_at(coordinate);
     }
-
-    Integer num_neighbors() const
-    {
-        return space.lock()->num_neighbors(coordinate);
-    }
-
-    Voxel get_neighbor(Integer nrand) const
-    {
-        return Voxel(world, space,
-                     space.lock()->get_neighbor(coordinate, nrand));
-    }
-
-    Voxel get_neighbor_randomly() const;
-
-    Voxel get_neighbor_randomly(Shape::dimension_kind dimension) const;
 };
 
 } // namespace spatiocyte
diff --git a/ecell4/spatiocyte/tests/SpatiocyteWorld_test.cpp b/ecell4/spatiocyte/tests/SpatiocyteWorld_test.cpp
index 369aa41cd..9df55759a 100644
--- a/ecell4/spatiocyte/tests/SpatiocyteWorld_test.cpp
+++ b/ecell4/spatiocyte/tests/SpatiocyteWorld_test.cpp
@@ -127,9 +127,9 @@ BOOST_AUTO_TEST_CASE(SpatiocyteWorld_test_neighbor)
     const Species sp("TEST", 1e-8, 1e-12);
     model->add_species_attribute(sp);
 
-    for (Integer i(0); i < voxel.num_neighbors(); ++i)
+    for (Integer i(0); i < world.num_neighbors(voxel); ++i)
     {
-        world.new_particle(sp, voxel.get_neighbor(i));
+        world.new_particle(sp, world.get_neighbor(voxel, i));
     }
     std::vector<std::pair<ParticleID, Particle>> particles(
         world.list_particles());

From 9f24487acd2d35ae9b7eb48b89571d3a53705c08 Mon Sep 17 00:00:00 2001
From: Suguru Kato <salmon.bioinformatics@gmail.com>
Date: Thu, 16 Jan 2020 13:26:52 +0900
Subject: [PATCH 06/21] :recycle: Change the types of intefaces_ and neighbors_

---
 ecell4/spatiocyte/SpatiocyteWorld.cpp | 37 ++++++++++++++++-----------
 ecell4/spatiocyte/SpatiocyteWorld.hpp |  4 +--
 ecell4/spatiocyte/Voxel.hpp           | 21 ++++++++++++++-
 3 files changed, 44 insertions(+), 18 deletions(-)

diff --git a/ecell4/spatiocyte/SpatiocyteWorld.cpp b/ecell4/spatiocyte/SpatiocyteWorld.cpp
index 56a423d95..f57d7b9c6 100644
--- a/ecell4/spatiocyte/SpatiocyteWorld.cpp
+++ b/ecell4/spatiocyte/SpatiocyteWorld.cpp
@@ -9,33 +9,40 @@ namespace ecell4
 namespace spatiocyte
 {
 
-void SpatiocyteWorld::add_space(std::unique_ptr<VoxelSpaceBase> space)
+void SpatiocyteWorld::add_space(std::unique_ptr<VoxelSpaceBase> uniq_space)
 {
+    boost::shared_ptr<VoxelSpaceBase> space(uniq_space.release());
+
     for (coordinate_type i(0); i < space->size(); ++i)
     {
-        const Real3 position(space->coordinate2position(i));
-        const coordinate_type nearest(
-            get_root()->position2coordinate(position));
+        const Voxel voxel(space, i);
+        const auto position(voxel.position());
+        const Voxel nearest(get_root(),
+                            get_root()->position2coordinate(position));
 
-        for (Integer j(0); j < get_root()->num_neighbors(nearest); ++j)
+        for (Integer j(0); j < num_neighbors(nearest); ++j)
         {
-            const coordinate_type neighbor(
-                get_root()->get_neighbor(nearest, j));
-            if (length(get_root()->coordinate2position(neighbor) - position) <
-                voxel_radius() * 2)
-                interfaces_.add(neighbor, i + size_);
+            const auto neighbor(get_neighbor(nearest, j));
+            if (length(neighbor.position() - position) < voxel_radius() * 2)
+            {
+                interfaces_.add(neighbor, voxel);
+            }
         }
     }
 
     for (const auto &interface : interfaces_)
     {
-        std::vector<coordinate_type> neighbors;
-        for (Integer i(0); i < get_root()->num_neighbors(interface.first); ++i)
+        const Voxel voxel(interface.first);
+
+        std::vector<Voxel> neighbors;
+        for (auto i(0); i < num_neighbors(voxel); ++i)
         {
-            const coordinate_type neighbor(
-                get_root()->get_neighbor(interface.first, i));
+            const auto neighbor(get_neighbor(voxel, i));
+
             if (!interfaces_.find(neighbor))
+            {
                 neighbors.push_back(neighbor);
+            }
         }
 
         for (const auto &adjoining : interface.second)
@@ -45,7 +52,7 @@ void SpatiocyteWorld::add_space(std::unique_ptr<VoxelSpaceBase> space)
     }
 
     size_ += space->size();
-    spaces_.push_back(space_type(space.release()));
+    spaces_.push_back(space);
 }
 
 void SpatiocyteWorld::set_value(const Species &sp, const Real value)
diff --git a/ecell4/spatiocyte/SpatiocyteWorld.hpp b/ecell4/spatiocyte/SpatiocyteWorld.hpp
index 25eac5432..5b00c5391 100644
--- a/ecell4/spatiocyte/SpatiocyteWorld.hpp
+++ b/ecell4/spatiocyte/SpatiocyteWorld.hpp
@@ -931,8 +931,8 @@ class SpatiocyteWorld : public WorldInterface
     std::size_t size_;
     space_container_type spaces_;
 
-    OneToManyMap<coordinate_type> interfaces_;
-    OneToManyMap<coordinate_type> neighbors_;
+    OneToManyMap<Voxel> interfaces_;
+    OneToManyMap<Voxel> neighbors_;
 
     boost::shared_ptr<RandomNumberGenerator> rng_;
     SerialIDGenerator<ParticleID> sidgen_;
diff --git a/ecell4/spatiocyte/Voxel.hpp b/ecell4/spatiocyte/Voxel.hpp
index 938997272..e4ad04602 100644
--- a/ecell4/spatiocyte/Voxel.hpp
+++ b/ecell4/spatiocyte/Voxel.hpp
@@ -3,6 +3,7 @@
 
 #include <boost/weak_ptr.hpp>
 #include <ecell4/core/VoxelSpaceBase.hpp>
+#include <ecell4/core/hash.hpp>
 #include <ecell4/core/types.hpp>
 
 namespace ecell4
@@ -15,7 +16,7 @@ class SpatiocyteWorld;
 
 struct Voxel
 {
-    typedef Integer coordinate_type;
+    typedef VoxelSpaceBase::coordinate_type coordinate_type;
 
     Voxel(boost::weak_ptr<VoxelSpaceBase> space, coordinate_type coordinate)
         : space(space), coordinate(coordinate)
@@ -37,10 +38,28 @@ struct Voxel
     {
         return space.lock()->get_voxel_pool_at(coordinate);
     }
+
+    bool operator==(const Voxel &rhs) const noexcept
+    {
+        return space.lock() == rhs.space.lock() && coordinate == rhs.coordinate;
+    }
 };
 
 } // namespace spatiocyte
 
 } // namespace ecell4
 
+ECELL4_DEFINE_HASH_BEGIN()
+template <>
+struct hash<ecell4::spatiocyte::Voxel>
+{
+    std::size_t operator()(const ecell4::spatiocyte::Voxel &val) const
+    {
+        auto ptr = val.space.lock().get();
+        return hash<decltype(ptr)>()(ptr) ^
+               static_cast<std::size_t>(val.coordinate);
+    }
+};
+ECELL4_DEFINE_HASH_END()
+
 #endif

From b4693e5e31f87b702d6a20bb3701c45c5e3182a7 Mon Sep 17 00:00:00 2001
From: Suguru Kato <salmon.bioinformatics@gmail.com>
Date: Thu, 16 Jan 2020 13:38:16 +0900
Subject: [PATCH 07/21] :recycle: Refactor add_structure() and remove
 is_inside()

---
 ecell4/spatiocyte/SpatiocyteWorld.cpp | 64 ++++++++++++++-------------
 ecell4/spatiocyte/SpatiocyteWorld.hpp |  5 ---
 2 files changed, 33 insertions(+), 36 deletions(-)

diff --git a/ecell4/spatiocyte/SpatiocyteWorld.cpp b/ecell4/spatiocyte/SpatiocyteWorld.cpp
index f57d7b9c6..2b98bb50b 100644
--- a/ecell4/spatiocyte/SpatiocyteWorld.cpp
+++ b/ecell4/spatiocyte/SpatiocyteWorld.cpp
@@ -234,26 +234,27 @@ SpatiocyteWorld::add_structure3(const Species &sp, const std::string &location,
                                 const boost::shared_ptr<const Shape> shape)
 {
     Integer count(0);
-    for (coordinate_type coord(0); coord < size(); ++coord)
+    for (const auto &space : spaces_)
     {
-        const Voxel voxel(coordinate2voxel(coord));
-
-        if (!this->is_inside(coord) || shape->is_inside(voxel.position()) > 0)
+        for (auto coord(0); coord < space->size(); ++coord)
         {
-            continue;
-        }
+            const Voxel voxel(space, coord);
+            if (shape->is_inside(voxel.position()) > 0)
+                continue;
 
-        if (voxel.get_voxel_pool()->species().serial() != location)
-        {
-            throw NotSupported("Mismatch in the location. Failed to place '" +
-                               sp.serial() + "' to '" +
-                               voxel.get_voxel_pool()->species().serial() +
-                               "'. " + "'" + location + "' is expected.");
-            continue;
-        }
+            if (voxel.get_voxel_pool()->species().serial() != location)
+            {
+                throw NotSupported(
+                    "Mismatch in the location. Failed to place '" +
+                    sp.serial() + "' to '" +
+                    voxel.get_voxel_pool()->species().serial() + "'. " + "'" +
+                    location + "' is expected.");
+                continue;
+            }
 
-        if (new_voxel_structure(sp, voxel))
-            ++count;
+            if (new_voxel_structure(sp, voxel))
+                ++count;
+        }
     }
     return count;
 }
@@ -263,25 +264,26 @@ SpatiocyteWorld::add_structure2(const Species &sp, const std::string &location,
                                 const boost::shared_ptr<const Shape> shape)
 {
     Integer count(0);
-    for (coordinate_type coord(0); coord < size(); ++coord)
+    for (const auto &space : spaces_)
     {
-        const Voxel voxel(coordinate2voxel(coord));
-        if (!this->is_inside(coord) || !is_surface_voxel(voxel, shape))
+        for (auto coord(0); coord < space->size(); ++coord)
         {
-            continue;
-        }
+            const Voxel voxel(space, coord);
+            if (!is_surface_voxel(voxel, shape))
+                continue;
+            if (voxel.get_voxel_pool()->species().serial() != location)
+            {
+                throw NotSupported(
+                    "Mismatch in the location. Failed to place '" +
+                    sp.serial() + "' to '" +
+                    voxel.get_voxel_pool()->species().serial() + "'. " + "'" +
+                    location + "' is expected.");
+                continue;
+            }
 
-        if (voxel.get_voxel_pool()->species().serial() != location)
-        {
-            throw NotSupported("Mismatch in the location. Failed to place '" +
-                               sp.serial() + "' to '" +
-                               voxel.get_voxel_pool()->species().serial() +
-                               "'. " + "'" + location + "' is expected.");
-            continue;
+            if (new_voxel_structure(sp, voxel))
+                ++count;
         }
-
-        if (new_voxel_structure(sp, voxel))
-            ++count;
     }
     return count;
 }
diff --git a/ecell4/spatiocyte/SpatiocyteWorld.hpp b/ecell4/spatiocyte/SpatiocyteWorld.hpp
index 5b00c5391..7f9d3accb 100644
--- a/ecell4/spatiocyte/SpatiocyteWorld.hpp
+++ b/ecell4/spatiocyte/SpatiocyteWorld.hpp
@@ -834,11 +834,6 @@ class SpatiocyteWorld : public WorldInterface
     }
 
 protected:
-    bool is_inside(const coordinate_type &coord) const
-    {
-        return get_root()->is_inside(coord);
-    }
-
     space_type get_root() const { return spaces_.at(0); }
 
     space_type get_space(coordinate_type &coordinate)

From 02cdd43e318adf358d0125db4338f8a030402456 Mon Sep 17 00:00:00 2001
From: Suguru Kato <salmon.bioinformatics@gmail.com>
Date: Thu, 16 Jan 2020 13:55:50 +0900
Subject: [PATCH 08/21] :recycle: Change the arguments of gen_particle_from()

---
 ecell4/spatiocyte/SpatiocyteWorld.hpp | 15 +++++++--------
 1 file changed, 7 insertions(+), 8 deletions(-)

diff --git a/ecell4/spatiocyte/SpatiocyteWorld.hpp b/ecell4/spatiocyte/SpatiocyteWorld.hpp
index 7f9d3accb..5d134b0ef 100644
--- a/ecell4/spatiocyte/SpatiocyteWorld.hpp
+++ b/ecell4/spatiocyte/SpatiocyteWorld.hpp
@@ -255,8 +255,7 @@ class SpatiocyteWorld : public WorldInterface
         {
             if (const auto &view = space->find_voxel(pid))
             {
-                return std::make_pair(
-                    pid, gen_particle_from(space, view->species, view->voxel));
+                return std::make_pair(pid, gen_particle_from(space, *view));
             }
         }
         throw "No particle corresponding to a given ParticleID is found.";
@@ -855,10 +854,12 @@ class SpatiocyteWorld : public WorldInterface
     bool is_surface_voxel(const Voxel &voxel,
                           const boost::shared_ptr<const Shape> shape) const;
 
-    Particle gen_particle_from(const space_type &space, const Species &species,
-                               const coordinate_type coordinate) const
+    Particle gen_particle_from(const space_type &space,
+                               const VoxelView &view) const
     {
-        const auto position = space->coordinate2position(coordinate);
+        const auto species = view.species;
+
+        const auto position = space->coordinate2position(view.voxel);
         const auto minfo_iter = molecule_info_cache_.find(species);
         if (minfo_iter != molecule_info_cache_.end())
         {
@@ -895,9 +896,7 @@ class SpatiocyteWorld : public WorldInterface
     {
         return map_voxels<std::pair<ParticleID, Particle>>(
             list_fn, [this](const space_type &space, const VoxelView &view) {
-                return std::make_pair(
-                    view.pid,
-                    gen_particle_from(space, view.species, view.voxel));
+                return std::make_pair(view.pid, gen_particle_from(space, view));
             });
     }
 

From cb60782a7df55702195f73f18d0798ac59979cb3 Mon Sep 17 00:00:00 2001
From: Suguru Kato <salmon.bioinformatics@gmail.com>
Date: Thu, 16 Jan 2020 14:13:41 +0900
Subject: [PATCH 09/21] :recycle: Add find_space_and_voxel_pool()

---
 ecell4/spatiocyte/ReactionEvent.cpp   | 11 +++++++----
 ecell4/spatiocyte/SpatiocyteWorld.hpp | 16 ++++++++++++++++
 2 files changed, 23 insertions(+), 4 deletions(-)

diff --git a/ecell4/spatiocyte/ReactionEvent.cpp b/ecell4/spatiocyte/ReactionEvent.cpp
index e89119934..eb84472b3 100644
--- a/ecell4/spatiocyte/ReactionEvent.cpp
+++ b/ecell4/spatiocyte/ReactionEvent.cpp
@@ -24,9 +24,12 @@ void ZerothOrderReactionEvent::fire_()
     {
         const MoleculeInfo info(world_->get_molecule_info(sp));
 
-        if (boost::shared_ptr<VoxelPool> location =
-                world_->find_voxel_pool(Species(info.loc)))
+        if (const auto space_and_location =
+                world_->find_space_and_voxel_pool(Species(info.loc)))
         {
+            const auto space = space_and_location->first;
+            const auto location = space_and_location->second;
+
             if (location->size() == 0)
             {
                 time_ += draw_dt();
@@ -35,8 +38,8 @@ void ZerothOrderReactionEvent::fire_()
 
             while (true)
             {
-                const Voxel voxel(world_->coordinate2voxel(
-                    world_->rng()->uniform_int(0, world_->size() - 1)));
+                const Voxel voxel(
+                    space, world_->rng()->uniform_int(0, world_->size() - 1));
 
                 if (voxel.get_voxel_pool() != location)
                 {
diff --git a/ecell4/spatiocyte/SpatiocyteWorld.hpp b/ecell4/spatiocyte/SpatiocyteWorld.hpp
index 5d134b0ef..dc135044a 100644
--- a/ecell4/spatiocyte/SpatiocyteWorld.hpp
+++ b/ecell4/spatiocyte/SpatiocyteWorld.hpp
@@ -407,6 +407,22 @@ class SpatiocyteWorld : public WorldInterface
         throw "No VoxelPool corresponding to a given Species is found";
     }
 
+    boost::optional<std::pair<space_type, boost::shared_ptr<const VoxelPool>>>
+    find_space_and_voxel_pool(const Species &species) const
+    {
+        for (const auto &space : spaces_)
+        {
+            if (space->has_species(species))
+            {
+                const auto voxel_pool = space->find_voxel_pool(species);
+                return std::pair<space_type,
+                                 boost::shared_ptr<const VoxelPool>>(
+                    space, voxel_pool);
+            }
+        }
+        return boost::none;
+    }
+
     bool has_molecule_pool(const Species &species) const
     {
         for (const auto &space : spaces_)

From 2e8516a5376367fae72136b3e09ea21b6b0722ed Mon Sep 17 00:00:00 2001
From: Suguru Kato <salmon.bioinformatics@gmail.com>
Date: Thu, 16 Jan 2020 14:39:08 +0900
Subject: [PATCH 10/21] :recycle: Add find_space_and_molecule_pool()

---
 ecell4/spatiocyte/SpatiocyteEvent.hpp | 23 ++++++++++++++---------
 ecell4/spatiocyte/SpatiocyteWorld.cpp | 26 +++++++++++++++-----------
 ecell4/spatiocyte/SpatiocyteWorld.hpp | 17 +++++++++++++++++
 3 files changed, 46 insertions(+), 20 deletions(-)

diff --git a/ecell4/spatiocyte/SpatiocyteEvent.hpp b/ecell4/spatiocyte/SpatiocyteEvent.hpp
index faf498c7a..a11bc2930 100644
--- a/ecell4/spatiocyte/SpatiocyteEvent.hpp
+++ b/ecell4/spatiocyte/SpatiocyteEvent.hpp
@@ -120,15 +120,20 @@ struct FirstOrderReactionEvent : SpatiocyteEvent
     ReactionInfo::Item choice()
     {
         const Species &species(rule_.reactants().at(0));
-        boost::shared_ptr<const MoleculePool> mt(
-            world_->find_molecule_pool(species));
-
-        const Integer i(rng_.lock()->uniform_int(0, mt->size() - 1));
-        const SpatiocyteWorld::coordinate_id_pair_type &info(mt->at(i));
-
-        // TODO: Calling coordinate2voxel() is invalid
-        return ReactionInfo::Item(info.pid, species,
-                                  world_->coordinate2voxel(info.coordinate));
+        if (const auto space_and_molecule_pool =
+                world_->find_space_and_molecule_pool(species))
+        {
+            const auto space = space_and_molecule_pool->first;
+            const auto molecule_pool = space_and_molecule_pool->second;
+
+            const auto i =
+                rng_.lock()->uniform_int(0, molecule_pool->size() - 1);
+            const auto &info = molecule_pool->at(i);
+
+            return ReactionInfo::Item(info.pid, species,
+                                      Voxel(space, info.coordinate));
+        }
+        throw "MoleculePool is not found";
     }
 
     boost::shared_ptr<SpatiocyteWorld> world_;
diff --git a/ecell4/spatiocyte/SpatiocyteWorld.cpp b/ecell4/spatiocyte/SpatiocyteWorld.cpp
index 2b98bb50b..33f0bcbf0 100644
--- a/ecell4/spatiocyte/SpatiocyteWorld.cpp
+++ b/ecell4/spatiocyte/SpatiocyteWorld.cpp
@@ -310,20 +310,24 @@ void SpatiocyteWorld::remove_molecules(const Species &sp, const Integer &num)
             "The number of molecules must be positive.");
     }
 
-    boost::shared_ptr<const MoleculePool> mtype(find_molecule_pool(sp));
-    if (mtype->size() < num)
+    if (const auto space_and_molecule_pool = find_space_and_molecule_pool(sp))
     {
-        throw std::invalid_argument(
-            "The number of molecules cannot be negative.");
-    }
+        const auto space = space_and_molecule_pool->first;
+        const auto molecule_pool = space_and_molecule_pool->second;
 
-    Integer count(0);
-    while (count < num)
-    {
-        const Integer idx(rng_->uniform_int(0, mtype->size() - 1));
-        if (coordinate2voxel(mtype->at(idx).coordinate).clear())
+        if (molecule_pool->size() < num)
+            throw std::invalid_argument(
+                "The number of molecules cannot be negative.");
+
+        auto count(0);
+        while (count < num)
         {
-            ++count;
+            const auto idx(rng_->uniform_int(0, molecule_pool->size() - 1));
+            const Voxel voxel(space, molecule_pool->at(idx).coordinate);
+            if (voxel.clear())
+            {
+                ++count;
+            }
         }
     }
 }
diff --git a/ecell4/spatiocyte/SpatiocyteWorld.hpp b/ecell4/spatiocyte/SpatiocyteWorld.hpp
index dc135044a..d28ceee70 100644
--- a/ecell4/spatiocyte/SpatiocyteWorld.hpp
+++ b/ecell4/spatiocyte/SpatiocyteWorld.hpp
@@ -454,6 +454,23 @@ class SpatiocyteWorld : public WorldInterface
         throw "No MoleculePool corresponding to a given Species is found";
     }
 
+    boost::optional<
+        std::pair<space_type, boost::shared_ptr<const MoleculePool>>>
+    find_space_and_molecule_pool(const Species &species) const
+    {
+        for (const auto &space : spaces_)
+        {
+            if (space->has_molecule_pool(species))
+            {
+                const auto molecule_pool(space->find_molecule_pool(species));
+                return std::pair<space_type,
+                                 boost::shared_ptr<const MoleculePool>>(
+                    space, molecule_pool);
+            }
+        }
+        return boost::none;
+    }
+
     /*
      * Coordinate Transformation
      */

From f6296178984dfad3ca5827a57614356ef9b36cab Mon Sep 17 00:00:00 2001
From: Suguru Kato <salmon.bioinformatics@gmail.com>
Date: Thu, 16 Jan 2020 14:54:48 +0900
Subject: [PATCH 11/21] :fire: Remove coordinate2voxel()

---
 ecell4/spatiocyte/ReactionEvent.cpp   |  2 +-
 ecell4/spatiocyte/SpatiocyteEvent.hpp |  1 +
 ecell4/spatiocyte/SpatiocyteWorld.cpp | 28 +++++++++++++++++----------
 ecell4/spatiocyte/SpatiocyteWorld.hpp | 19 ++++--------------
 ecell4/spatiocyte/StepEvent.cpp       | 22 ++++++++++++++-------
 5 files changed, 39 insertions(+), 33 deletions(-)

diff --git a/ecell4/spatiocyte/ReactionEvent.cpp b/ecell4/spatiocyte/ReactionEvent.cpp
index eb84472b3..4836fcb5c 100644
--- a/ecell4/spatiocyte/ReactionEvent.cpp
+++ b/ecell4/spatiocyte/ReactionEvent.cpp
@@ -39,7 +39,7 @@ void ZerothOrderReactionEvent::fire_()
             while (true)
             {
                 const Voxel voxel(
-                    space, world_->rng()->uniform_int(0, world_->size() - 1));
+                    space, world_->rng()->uniform_int(0, space->size() - 1));
 
                 if (voxel.get_voxel_pool() != location)
                 {
diff --git a/ecell4/spatiocyte/SpatiocyteEvent.hpp b/ecell4/spatiocyte/SpatiocyteEvent.hpp
index a11bc2930..5ecd0e0ee 100644
--- a/ecell4/spatiocyte/SpatiocyteEvent.hpp
+++ b/ecell4/spatiocyte/SpatiocyteEvent.hpp
@@ -66,6 +66,7 @@ struct StepEvent : SpatiocyteEvent
 protected:
     boost::shared_ptr<Model> model_;
     boost::shared_ptr<SpatiocyteWorld> world_;
+    boost::weak_ptr<VoxelSpaceBase> space_;
     boost::shared_ptr<MoleculePool> mpool_;
 
     const Real alpha_;
diff --git a/ecell4/spatiocyte/SpatiocyteWorld.cpp b/ecell4/spatiocyte/SpatiocyteWorld.cpp
index 33f0bcbf0..f91dc29ae 100644
--- a/ecell4/spatiocyte/SpatiocyteWorld.cpp
+++ b/ecell4/spatiocyte/SpatiocyteWorld.cpp
@@ -156,21 +156,29 @@ bool SpatiocyteWorld::add_molecules(const Species &sp, const Integer &num)
 
     const MoleculeInfo info(get_molecule_info(sp));
 
-    Integer count(0);
-    while (count < num)
+    if (const auto space_and_location =
+            find_space_and_voxel_pool(Species(info.loc)))
     {
-        const Voxel voxel(coordinate2voxel(rng()->uniform_int(0, size() - 1)));
+        const auto space = space_and_location->first;
+        const auto location = space_and_location->second;
 
-        if (voxel.get_voxel_pool()->species().serial() != info.loc)
-        {
-            continue;
-        }
-        else if (new_particle(sp, voxel))
+        if (location->size() < num)
+            return false;
+
+        auto count(0);
+        while (count < num)
         {
-            ++count;
+            const Voxel voxel(space, rng()->uniform_int(0, space->size() - 1));
+
+            if (voxel.get_voxel_pool()->species().serial() != info.loc)
+                continue;
+
+            if (new_particle(sp, voxel))
+                ++count;
         }
+        return true;
     }
-    return true;
+    return false;
 }
 
 bool SpatiocyteWorld::add_molecules(const Species &sp, const Integer &num,
diff --git a/ecell4/spatiocyte/SpatiocyteWorld.hpp b/ecell4/spatiocyte/SpatiocyteWorld.hpp
index d28ceee70..ce9fa04c1 100644
--- a/ecell4/spatiocyte/SpatiocyteWorld.hpp
+++ b/ecell4/spatiocyte/SpatiocyteWorld.hpp
@@ -407,7 +407,7 @@ class SpatiocyteWorld : public WorldInterface
         throw "No VoxelPool corresponding to a given Species is found";
     }
 
-    boost::optional<std::pair<space_type, boost::shared_ptr<const VoxelPool>>>
+    boost::optional<std::pair<space_type, boost::shared_ptr<VoxelPool>>>
     find_space_and_voxel_pool(const Species &species) const
     {
         for (const auto &space : spaces_)
@@ -415,8 +415,7 @@ class SpatiocyteWorld : public WorldInterface
             if (space->has_species(species))
             {
                 const auto voxel_pool = space->find_voxel_pool(species);
-                return std::pair<space_type,
-                                 boost::shared_ptr<const VoxelPool>>(
+                return std::pair<space_type, boost::shared_ptr<VoxelPool>>(
                     space, voxel_pool);
             }
         }
@@ -454,8 +453,7 @@ class SpatiocyteWorld : public WorldInterface
         throw "No MoleculePool corresponding to a given Species is found";
     }
 
-    boost::optional<
-        std::pair<space_type, boost::shared_ptr<const MoleculePool>>>
+    boost::optional<std::pair<space_type, boost::shared_ptr<MoleculePool>>>
     find_space_and_molecule_pool(const Species &species) const
     {
         for (const auto &space : spaces_)
@@ -463,8 +461,7 @@ class SpatiocyteWorld : public WorldInterface
             if (space->has_molecule_pool(species))
             {
                 const auto molecule_pool(space->find_molecule_pool(species));
-                return std::pair<space_type,
-                                 boost::shared_ptr<const MoleculePool>>(
+                return std::pair<space_type, boost::shared_ptr<MoleculePool>>(
                     space, molecule_pool);
             }
         }
@@ -943,14 +940,6 @@ class SpatiocyteWorld : public WorldInterface
             });
     }
 
-public:
-    // TODO: Calling this function is invalid, and this should be removed.
-    Voxel coordinate2voxel(const coordinate_type &coordinate)
-    {
-        coordinate_type coord(coordinate);
-        return Voxel(get_space(coord), coord);
-    }
-
 protected:
     typedef utils::get_mapper_mf<Species, MoleculeInfo>::type
         molecule_info_cache_t;
diff --git a/ecell4/spatiocyte/StepEvent.cpp b/ecell4/spatiocyte/StepEvent.cpp
index f5cad1e04..65e0c2fd0 100644
--- a/ecell4/spatiocyte/StepEvent.cpp
+++ b/ecell4/spatiocyte/StepEvent.cpp
@@ -10,9 +10,19 @@ namespace spatiocyte
 StepEvent::StepEvent(boost::shared_ptr<Model> model,
                      boost::shared_ptr<SpatiocyteWorld> world,
                      const Species &species, const Real &t, const Real alpha)
-    : SpatiocyteEvent(t), model_(model), world_(world),
-      mpool_(world_->find_molecule_pool(species)), alpha_(alpha)
+    : SpatiocyteEvent(t), model_(model), world_(world), alpha_(alpha)
 {
+    if (const auto space_and_molecule_pool =
+            world_->find_space_and_molecule_pool(species))
+    {
+        space_ = space_and_molecule_pool->first;
+        mpool_ = space_and_molecule_pool->second;
+    }
+    else
+    {
+        throw "MoleculePool is not found";
+    }
+
     time_ = t;
 }
 
@@ -48,7 +58,7 @@ void StepEvent3D::walk(const Real &alpha)
     std::size_t idx(0);
     for (const auto &info : voxels)
     {
-        const Voxel voxel(world_->coordinate2voxel(info.coordinate));
+        const Voxel voxel(space_, info.coordinate);
 
         if (voxel.get_voxel_pool() != mpool_)
         {
@@ -105,8 +115,7 @@ void StepEvent2D::walk(const Real &alpha)
     std::size_t idx(0);
     for (const auto &info : voxels)
     {
-        // TODO: Calling coordinate2voxel is invalid
-        const Voxel voxel(world_->coordinate2voxel(info.coordinate));
+        const Voxel voxel(space_, info.coordinate);
 
         if (voxel.get_voxel_pool() != mpool_)
         {
@@ -135,8 +144,7 @@ void StepEvent::attempt_reaction_(
     const SpatiocyteWorld::coordinate_id_pair_type &info, const Voxel &dst,
     const Real &alpha)
 {
-    // TODO: Calling coordiante2voxel is invalid
-    const Voxel voxel(world_->coordinate2voxel(info.coordinate));
+    const Voxel voxel(space_, info.coordinate);
     boost::shared_ptr<const VoxelPool> from_mt(voxel.get_voxel_pool());
     boost::shared_ptr<const VoxelPool> to_mt(dst.get_voxel_pool());
 

From b889161d5301741d77f204669a623c38b93c2f3b Mon Sep 17 00:00:00 2001
From: Suguru Kato <salmon.bioinformatics@gmail.com>
Date: Thu, 16 Jan 2020 15:56:53 +0900
Subject: [PATCH 12/21] :bug: Fix bug around locating molecules

---
 ecell4/core/VoxelSpaceBase.cpp        | 6 ++++++
 ecell4/spatiocyte/SpatiocyteWorld.cpp | 7 +++++--
 2 files changed, 11 insertions(+), 2 deletions(-)

diff --git a/ecell4/core/VoxelSpaceBase.cpp b/ecell4/core/VoxelSpaceBase.cpp
index 941d4b55a..c29d330cd 100644
--- a/ecell4/core/VoxelSpaceBase.cpp
+++ b/ecell4/core/VoxelSpaceBase.cpp
@@ -62,6 +62,12 @@ bool VoxelSpaceBase::has_voxel(const ParticleID &pid) const
 
 Integer VoxelSpaceBase::num_voxels_exact(const Species &sp) const
 {
+    // TODO: change the way of dealing vacant
+    if (sp.serial() == "")
+    {
+        return vacant_->size();
+    }
+
     {
         voxel_pool_map_type::const_iterator itr(voxel_pools_.find(sp));
         if (itr != voxel_pools_.end())
diff --git a/ecell4/spatiocyte/SpatiocyteWorld.cpp b/ecell4/spatiocyte/SpatiocyteWorld.cpp
index f91dc29ae..091a49ee8 100644
--- a/ecell4/spatiocyte/SpatiocyteWorld.cpp
+++ b/ecell4/spatiocyte/SpatiocyteWorld.cpp
@@ -247,7 +247,9 @@ SpatiocyteWorld::add_structure3(const Species &sp, const std::string &location,
         for (auto coord(0); coord < space->size(); ++coord)
         {
             const Voxel voxel(space, coord);
-            if (shape->is_inside(voxel.position()) > 0)
+            // should check if coord doesn't point at the boundary.
+            if (!space->is_inside(coord) ||
+                shape->is_inside(voxel.position()) > 0)
                 continue;
 
             if (voxel.get_voxel_pool()->species().serial() != location)
@@ -277,7 +279,8 @@ SpatiocyteWorld::add_structure2(const Species &sp, const std::string &location,
         for (auto coord(0); coord < space->size(); ++coord)
         {
             const Voxel voxel(space, coord);
-            if (!is_surface_voxel(voxel, shape))
+            // should check if coord doesn't point at the boundary.
+            if (!space->is_inside(coord) || !is_surface_voxel(voxel, shape))
                 continue;
             if (voxel.get_voxel_pool()->species().serial() != location)
             {

From 4759c2de00db51303b771820f3b2e0946a6150c2 Mon Sep 17 00:00:00 2001
From: Suguru Kato <salmon.bioinformatics@gmail.com>
Date: Thu, 16 Jan 2020 16:00:40 +0900
Subject: [PATCH 13/21] :fire: Remove SpatiocyteWorld::coordinate_type

---
 ecell4/spatiocyte/SpatiocyteWorld.cpp |  2 +-
 ecell4/spatiocyte/SpatiocyteWorld.hpp | 16 ----------------
 2 files changed, 1 insertion(+), 17 deletions(-)

diff --git a/ecell4/spatiocyte/SpatiocyteWorld.cpp b/ecell4/spatiocyte/SpatiocyteWorld.cpp
index 091a49ee8..770305a46 100644
--- a/ecell4/spatiocyte/SpatiocyteWorld.cpp
+++ b/ecell4/spatiocyte/SpatiocyteWorld.cpp
@@ -13,7 +13,7 @@ void SpatiocyteWorld::add_space(std::unique_ptr<VoxelSpaceBase> uniq_space)
 {
     boost::shared_ptr<VoxelSpaceBase> space(uniq_space.release());
 
-    for (coordinate_type i(0); i < space->size(); ++i)
+    for (auto i(0); i < space->size(); ++i)
     {
         const Voxel voxel(space, i);
         const auto position(voxel.position());
diff --git a/ecell4/spatiocyte/SpatiocyteWorld.hpp b/ecell4/spatiocyte/SpatiocyteWorld.hpp
index ce9fa04c1..b655387a3 100644
--- a/ecell4/spatiocyte/SpatiocyteWorld.hpp
+++ b/ecell4/spatiocyte/SpatiocyteWorld.hpp
@@ -42,7 +42,6 @@ class SpatiocyteWorld : public WorldInterface
     typedef LatticeSpaceVectorImpl default_root_type;
 
     typedef VoxelSpaceBase::coordinate_id_pair_type coordinate_id_pair_type;
-    typedef VoxelSpaceBase::coordinate_type coordinate_type;
 
     typedef boost::shared_ptr<VoxelSpaceBase> space_type;
     typedef std::vector<space_type> space_container_type;
@@ -725,7 +724,6 @@ class SpatiocyteWorld : public WorldInterface
 
     std::vector<Species> list_non_structure_species() const;
     std::vector<Species> list_structure_species() const;
-    // std::vector<coordinate_type> list_coords(const Species& sp) const;
 
     boost::optional<ParticleID> new_particle(const Species &sp,
                                              const Voxel &voxel)
@@ -790,8 +788,6 @@ class SpatiocyteWorld : public WorldInterface
     const Voxel get_neighbor_randomly(const Voxel &voxel,
                                       Shape::dimension_kind dimension) const;
 
-    // bool update_molecule(coordinate_type at, Species species);
-
     const Species &draw_species(const Species &pttrn) const;
 
     boost::shared_ptr<RandomNumberGenerator> rng() const { return rng_; }
@@ -865,18 +861,6 @@ class SpatiocyteWorld : public WorldInterface
 protected:
     space_type get_root() const { return spaces_.at(0); }
 
-    space_type get_space(coordinate_type &coordinate)
-    {
-        for (const auto &space : spaces_)
-        {
-            if (coordinate < space->size())
-                return space;
-
-            coordinate -= space->size();
-        }
-        return space_type();
-    }
-
     Integer add_structure2(const Species &sp, const std::string &location,
                            const boost::shared_ptr<const Shape> shape);
     Integer add_structure3(const Species &sp, const std::string &location,

From d99f989a083dd2d6e2dcd6b8a383684f5dd74da1 Mon Sep 17 00:00:00 2001
From: Suguru Kato <salmon.bioinformatics@gmail.com>
Date: Thu, 23 Jan 2020 11:07:37 +0900
Subject: [PATCH 14/21] :bug: Fix num_voxels* and num_molecules*

---
 ecell4/core/VoxelSpaceBase.cpp             | 21 +++++++++++++++++++--
 ecell4/core/tests/LatticeSpace_test.cpp    |  2 ++
 ecell4/core/tests/OffLatticeSpace_test.cpp |  2 +-
 3 files changed, 22 insertions(+), 3 deletions(-)

diff --git a/ecell4/core/VoxelSpaceBase.cpp b/ecell4/core/VoxelSpaceBase.cpp
index c29d330cd..559a2f8d2 100644
--- a/ecell4/core/VoxelSpaceBase.cpp
+++ b/ecell4/core/VoxelSpaceBase.cpp
@@ -21,6 +21,14 @@ Integer VoxelSpaceBase::num_molecules(const Species &sp) const
     Integer count(0);
     SpeciesExpressionMatcher sexp(sp);
 
+    {
+        const auto cnt = sexp.count(vacant_->species());
+        if (cnt > 0)
+        {
+            count += vacant_->size() * cnt;
+        }
+    }
+
     for (voxel_pool_map_type::const_iterator itr(voxel_pools_.begin());
          itr != voxel_pools_.end(); ++itr)
     {
@@ -62,8 +70,7 @@ bool VoxelSpaceBase::has_voxel(const ParticleID &pid) const
 
 Integer VoxelSpaceBase::num_voxels_exact(const Species &sp) const
 {
-    // TODO: change the way of dealing vacant
-    if (sp.serial() == "")
+    if (sp == vacant_->species())
     {
         return vacant_->size();
     }
@@ -93,6 +100,11 @@ Integer VoxelSpaceBase::num_voxels(const Species &sp) const
     Integer count(0);
     SpeciesExpressionMatcher sexp(sp);
 
+    if (sexp.match(vacant_->species()))
+    {
+        count += vacant_->size();
+    }
+
     for (voxel_pool_map_type::const_iterator itr(voxel_pools_.begin());
          itr != voxel_pools_.end(); ++itr)
     {
@@ -118,6 +130,11 @@ Integer VoxelSpaceBase::num_voxels() const
 {
     Integer count(0);
 
+    if (vacant_->species().serial() != "")
+    {
+        count += vacant_->size();
+    }
+
     for (voxel_pool_map_type::const_iterator itr(voxel_pools_.begin());
          itr != voxel_pools_.end(); ++itr)
     {
diff --git a/ecell4/core/tests/LatticeSpace_test.cpp b/ecell4/core/tests/LatticeSpace_test.cpp
index 0b60b11d2..c76855f7c 100644
--- a/ecell4/core/tests/LatticeSpace_test.cpp
+++ b/ecell4/core/tests/LatticeSpace_test.cpp
@@ -40,6 +40,8 @@ BOOST_AUTO_TEST_CASE(LatticeSpace_test_constructor) { ; }
 BOOST_AUTO_TEST_CASE(CheckVacantSize)
 {
     BOOST_CHECK_EQUAL(space.actual_size(), space.vacant()->size());
+    BOOST_CHECK_EQUAL(space.num_voxels_exact(Species("")),
+                      space.vacant()->size());
 }
 
 BOOST_AUTO_TEST_CASE(GetVoxel)
diff --git a/ecell4/core/tests/OffLatticeSpace_test.cpp b/ecell4/core/tests/OffLatticeSpace_test.cpp
index 2746355fb..a440813db 100644
--- a/ecell4/core/tests/OffLatticeSpace_test.cpp
+++ b/ecell4/core/tests/OffLatticeSpace_test.cpp
@@ -70,7 +70,7 @@ BOOST_AUTO_TEST_CASE(OffLatticeSpace_test_voxelspacebase)
     BOOST_CHECK_EQUAL(space.list_species().size(), 1);
     BOOST_CHECK_EQUAL(space.num_voxels_exact(species), 1);
     BOOST_CHECK_EQUAL(space.num_voxels(species), 1);
-    BOOST_CHECK_EQUAL(space.num_voxels(), 1);
+    BOOST_CHECK_EQUAL(space.num_voxels(), 10);
 
     BOOST_CHECK(space.has_voxel(pid));
     BOOST_CHECK(!space.has_voxel(sidgen()));

From de79eb5dab627e40054a79fa93ccc6ceda5024f7 Mon Sep 17 00:00:00 2001
From: Suguru Kato <salmon.bioinformatics@gmail.com>
Date: Thu, 23 Jan 2020 12:03:05 +0900
Subject: [PATCH 15/21] :bug: Fix the bug of SpatiocyteWorld::add_molecules()

---
 ecell4/core/VoxelSpaceBase.cpp             |  4 ++--
 ecell4/core/tests/OffLatticeSpace_test.cpp |  2 +-
 ecell4/spatiocyte/SpatiocyteWorld.cpp      |  2 +-
 tests/spatiocyte/offlattice.py             | 15 +++++++++++++++
 4 files changed, 19 insertions(+), 4 deletions(-)

diff --git a/ecell4/core/VoxelSpaceBase.cpp b/ecell4/core/VoxelSpaceBase.cpp
index 559a2f8d2..a799e3b64 100644
--- a/ecell4/core/VoxelSpaceBase.cpp
+++ b/ecell4/core/VoxelSpaceBase.cpp
@@ -202,7 +202,7 @@ VoxelSpaceBase::list_voxels_exact(const Species &sp) const
 
 boost::shared_ptr<VoxelPool> VoxelSpaceBase::find_voxel_pool(const Species &sp)
 {
-    if (sp.serial() == "")
+    if (sp == vacant_->species())
         return vacant_;
 
     voxel_pool_map_type::iterator itr(voxel_pools_.find(sp));
@@ -216,7 +216,7 @@ boost::shared_ptr<VoxelPool> VoxelSpaceBase::find_voxel_pool(const Species &sp)
 boost::shared_ptr<const VoxelPool>
 VoxelSpaceBase::find_voxel_pool(const Species &sp) const
 {
-    if (sp.serial() == "")
+    if (sp == vacant_->species())
         return vacant_;
 
     voxel_pool_map_type::const_iterator itr(voxel_pools_.find(sp));
diff --git a/ecell4/core/tests/OffLatticeSpace_test.cpp b/ecell4/core/tests/OffLatticeSpace_test.cpp
index a440813db..d580ce7f1 100644
--- a/ecell4/core/tests/OffLatticeSpace_test.cpp
+++ b/ecell4/core/tests/OffLatticeSpace_test.cpp
@@ -41,7 +41,7 @@ struct Fixture
         for (int i(1); i < 10; ++i)
             adjoining_pairs.push_back(std::make_pair(i - 1, i));
         space = OffLatticeSpace(voxel_radius, base, positions, adjoining_pairs);
-        space.make_molecular_type(species, "");
+        space.make_molecular_type(species, "Base");
     }
 };
 
diff --git a/ecell4/spatiocyte/SpatiocyteWorld.cpp b/ecell4/spatiocyte/SpatiocyteWorld.cpp
index 770305a46..c62f2dbaf 100644
--- a/ecell4/spatiocyte/SpatiocyteWorld.cpp
+++ b/ecell4/spatiocyte/SpatiocyteWorld.cpp
@@ -170,7 +170,7 @@ bool SpatiocyteWorld::add_molecules(const Species &sp, const Integer &num)
         {
             const Voxel voxel(space, rng()->uniform_int(0, space->size() - 1));
 
-            if (voxel.get_voxel_pool()->species().serial() != info.loc)
+            if (voxel.get_voxel_pool() != location)
                 continue;
 
             if (new_particle(sp, voxel))
diff --git a/tests/spatiocyte/offlattice.py b/tests/spatiocyte/offlattice.py
index ccfa9d8b9..105aed8c1 100644
--- a/tests/spatiocyte/offlattice.py
+++ b/tests/spatiocyte/offlattice.py
@@ -15,3 +15,18 @@ def test_constructor(self):
         species = Species('Base')
         world = SpatiocyteWorld(ones(), self.voxel_radius)
         world.add_offlattice(species, self.offlattice)
+
+    def test_add_molecules(self):
+        base = Species('Base', radius=self.voxel_radius, D=0.0)
+        species = Species('A', radius=self.voxel_radius, D=1.0, location='Base')
+
+        model = NetworkModel()
+        model.add_species_attribute(base)
+        model.add_species_attribute(species)
+
+        world = SpatiocyteWorld(ones(), self.voxel_radius)
+        world.bind_to(model)
+        world.add_offlattice(base, self.offlattice)
+
+        world.add_molecules(species, 5)
+        self.assertEqual(5, world.num_molecules(species))

From 12ac0c81bd4401a80080af1e1b45c94baa43527d Mon Sep 17 00:00:00 2001
From: Suguru Kato <salmon.bioinformatics@gmail.com>
Date: Thu, 23 Jan 2020 13:36:17 +0900
Subject: [PATCH 16/21] :recycle: Refactor get_molecule_info()

---
 ecell4/spatiocyte/SpatiocyteWorld.hpp | 65 ++++++++++-----------------
 1 file changed, 24 insertions(+), 41 deletions(-)

diff --git a/ecell4/spatiocyte/SpatiocyteWorld.hpp b/ecell4/spatiocyte/SpatiocyteWorld.hpp
index b655387a3..52c10c779 100644
--- a/ecell4/spatiocyte/SpatiocyteWorld.hpp
+++ b/ecell4/spatiocyte/SpatiocyteWorld.hpp
@@ -30,10 +30,10 @@ namespace spatiocyte
 
 struct MoleculeInfo
 {
-    const Real radius;
-    const Real D;
-    const std::string loc;
-    const Shape::dimension_kind dimension;
+    Real radius;
+    Real D;
+    std::string loc;
+    Shape::dimension_kind dimension;
 };
 
 class SpatiocyteWorld : public WorldInterface
@@ -559,11 +559,13 @@ class SpatiocyteWorld : public WorldInterface
      */
     const MoleculeInfo get_molecule_info(const Species &sp)
     {
-        boost::optional<Real> diffusion_coef;
-        boost::optional<Real> radius;
-        boost::optional<std::string> location;
-        Shape::dimension_kind dimension =
-            Shape::THREE; // The default dimension is 3
+        // Default
+        MoleculeInfo info = {
+            /* radius = */ voxel_radius(),
+            /* D = */ 0.0,
+            /* loc = */ "",
+            /* dimension = */ Shape::THREE,
+        };
 
         const auto itr = molecule_info_cache_.find(sp);
         if (itr != molecule_info_cache_.end())
@@ -571,10 +573,7 @@ class SpatiocyteWorld : public WorldInterface
             // return itr->second;
             // TODO: the below code is only for warning.
             //       In the future, the value should be returned immediately.
-            diffusion_coef = itr->second.D;
-            radius = itr->second.radius;
-            location = itr->second.loc;
-            dimension = itr->second.dimension;
+            info = itr->second;
         }
         else if (const auto model = lock_model())
         {
@@ -582,70 +581,54 @@ class SpatiocyteWorld : public WorldInterface
 
             if (species_from_model.has_attribute("D"))
             {
-                diffusion_coef = species_from_model.get_attribute_as<Real>("D");
+                info.D = species_from_model.get_attribute_as<Real>("D");
             }
 
             if (species_from_model.has_attribute("radius"))
             {
-                radius = species_from_model.get_attribute_as<Real>("radius");
+                info.radius =
+                    species_from_model.get_attribute_as<Real>("radius");
             }
 
             if (species_from_model.has_attribute("location"))
             {
-                location = species_from_model.get_attribute_as<std::string>(
+                info.loc = species_from_model.get_attribute_as<std::string>(
                     "location");
             }
 
-            dimension = extras::get_dimension_from_model(sp, model);
+            info.dimension = extras::get_dimension_from_model(sp, model);
         }
 
         if (sp.has_attribute("D"))
         {
             const auto new_value = sp.get_attribute_as<Real>("D");
-            if (diffusion_coef && *diffusion_coef != new_value)
+            if (info.D != new_value)
             {
                 warning("D");
-                diffusion_coef = new_value;
+                info.D = new_value;
             }
         }
 
         if (sp.has_attribute("radius"))
         {
             const auto new_value = sp.get_attribute_as<Real>("radius");
-            if (radius && *radius != new_value)
+            if (info.radius != new_value)
             {
                 warning("radius");
-                radius = new_value;
+                info.radius = new_value;
             }
         }
 
         if (sp.has_attribute("location"))
         {
             const auto new_value = sp.get_attribute_as<std::string>("location");
-            if (location && *location != new_value)
+            if (info.loc != new_value)
             {
                 warning("location");
-                location = new_value;
+                info.loc = new_value;
             }
         }
 
-        if (diffusion_coef == boost::none)
-        {
-            diffusion_coef = 0.0;
-        }
-
-        if (radius == boost::none)
-        {
-            radius = voxel_radius();
-        }
-
-        if (location == boost::none)
-        {
-            location = "";
-        }
-
-        const MoleculeInfo info = {*radius, *diffusion_coef, *location,
-                                   dimension};
         molecule_info_cache_.insert(
             molecule_info_cache_t::value_type(sp, info));
         return info;
@@ -939,7 +922,7 @@ class SpatiocyteWorld : public WorldInterface
 
     boost::weak_ptr<Model> model_;
     molecule_info_cache_t molecule_info_cache_;
-};
+}; // namespace spatiocyte
 
 inline SpatiocyteWorld *create_spatiocyte_world_cell_list_impl(
     const Real3 &edge_lengths, const Real &voxel_radius,

From aef1f5f5b4f361f504dc9683e74b466b46d4ed1b Mon Sep 17 00:00:00 2001
From: Suguru Kato <salmon.bioinformatics@gmail.com>
Date: Thu, 23 Jan 2020 13:36:57 +0900
Subject: [PATCH 17/21] :bug: Get information from a model before add a space

---
 ecell4/python_api/spatiocyte.cpp | 5 ++++-
 1 file changed, 4 insertions(+), 1 deletion(-)

diff --git a/ecell4/python_api/spatiocyte.cpp b/ecell4/python_api/spatiocyte.cpp
index c73318584..4adc77b50 100644
--- a/ecell4/python_api/spatiocyte.cpp
+++ b/ecell4/python_api/spatiocyte.cpp
@@ -155,7 +155,10 @@ static inline void define_spatiocyte_world(py::module &m)
         .def("add_offlattice",
              [](SpatiocyteWorld &self, const Species &species,
                 const OffLattice &offlattice) {
-                 self.add_space(offlattice.generate_space(species));
+                 const auto info = self.get_molecule_info(species);
+                 const auto updated = Species(species.serial(), info.radius,
+                                              info.D, info.loc, info.dimension);
+                 self.add_space(offlattice.generate_space(updated));
              })
         .def_static("calculate_voxel_volume",
                     &SpatiocyteWorld::calculate_voxel_volume)

From e35613160384e1dfe7946b5e077a3873a4c347d2 Mon Sep 17 00:00:00 2001
From: Suguru Kato <salmon.bioinformatics@gmail.com>
Date: Thu, 23 Jan 2020 13:38:13 +0900
Subject: [PATCH 18/21] :white_check_mark: Add test for get_molecule_info()

---
 .../spatiocyte/tests/SpatiocyteWorld_test.cpp | 22 ++++++++++++++++++-
 1 file changed, 21 insertions(+), 1 deletion(-)

diff --git a/ecell4/spatiocyte/tests/SpatiocyteWorld_test.cpp b/ecell4/spatiocyte/tests/SpatiocyteWorld_test.cpp
index 9df55759a..6a8c433b6 100644
--- a/ecell4/spatiocyte/tests/SpatiocyteWorld_test.cpp
+++ b/ecell4/spatiocyte/tests/SpatiocyteWorld_test.cpp
@@ -27,7 +27,7 @@ struct Fixture
 
     Fixture()
         : edge_lengths(1e-6, 1e-6, 1e-6), voxel_radius(1e-8),
-          rng(new GSLRandomNumberGenerator()), model(new NetworkModel),
+          rng(new GSLRandomNumberGenerator()), model(new NetworkModel()),
           world(edge_lengths, voxel_radius, rng)
     {
         world.bind_to(model);
@@ -63,6 +63,26 @@ BOOST_AUTO_TEST_CASE(SpatiocyteWorld_test_list_particles)
     BOOST_CHECK_EQUAL(particles.size(), 0);
 }
 
+BOOST_AUTO_TEST_CASE(SpatiocyteWorld_test_get_molecule_info)
+{
+    const Species m("M", voxel_radius, 0.0);
+    const Species a("A", voxel_radius, 1.0, "M");
+    model->add_species_attribute(m);
+    model->add_species_attribute(a);
+
+    const auto info_m = world.get_molecule_info(m);
+    BOOST_CHECK_EQUAL(info_m.radius, voxel_radius);
+    BOOST_CHECK_EQUAL(info_m.D, 0.0);
+    BOOST_CHECK_EQUAL(info_m.loc, "");
+    BOOST_CHECK_EQUAL(info_m.dimension, Shape::THREE);
+
+    const auto info_a = world.get_molecule_info(a);
+    BOOST_CHECK_EQUAL(info_a.radius, voxel_radius);
+    BOOST_CHECK_EQUAL(info_a.D, 1.0);
+    BOOST_CHECK_EQUAL(info_a.loc, "M");
+    BOOST_CHECK_EQUAL(info_a.dimension, Shape::THREE);
+}
+
 BOOST_AUTO_TEST_CASE(SpatiocyteWorld_test_update_particles)
 {
     SerialIDGenerator<ParticleID> sidgen;

From 186b1ecbdfbf8d26728fb6821d08fa9a3faad777 Mon Sep 17 00:00:00 2001
From: Suguru Kato <salmon.bioinformatics@gmail.com>
Date: Thu, 23 Jan 2020 13:39:13 +0900
Subject: [PATCH 19/21] :bug: Fix get_location() with a vacant voxel

---
 ecell4/spatiocyte/SpatiocyteReactions.cpp | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/ecell4/spatiocyte/SpatiocyteReactions.cpp b/ecell4/spatiocyte/SpatiocyteReactions.cpp
index 472c55388..95dc3d7f8 100644
--- a/ecell4/spatiocyte/SpatiocyteReactions.cpp
+++ b/ecell4/spatiocyte/SpatiocyteReactions.cpp
@@ -23,7 +23,7 @@ inline const std::string get_location(boost::shared_ptr<SpatiocyteWorld> world,
     if (mtype->is_vacant())
         return "";
     boost::shared_ptr<const VoxelPool> ltype(mtype->location());
-    return ltype->is_vacant() ? "" : ltype->species().serial();
+    return ltype->species().serial();
 }
 
 static inline void make_product(boost::shared_ptr<SpatiocyteWorld> world,

From b0fe1aaa172f960a0460f2fde9a357b08d98eff7 Mon Sep 17 00:00:00 2001
From: Suguru Kato <salmon.bioinformatics@gmail.com>
Date: Thu, 23 Jan 2020 14:00:26 +0900
Subject: [PATCH 20/21] :bug: Fix SpatiocyteWorld::check_neighbor()

---
 ecell4/spatiocyte/SpatiocyteWorld.cpp | 3 +--
 1 file changed, 1 insertion(+), 2 deletions(-)

diff --git a/ecell4/spatiocyte/SpatiocyteWorld.cpp b/ecell4/spatiocyte/SpatiocyteWorld.cpp
index c62f2dbaf..e286be8f7 100644
--- a/ecell4/spatiocyte/SpatiocyteWorld.cpp
+++ b/ecell4/spatiocyte/SpatiocyteWorld.cpp
@@ -355,8 +355,7 @@ boost::optional<Voxel> SpatiocyteWorld::check_neighbor(const Voxel &voxel,
     {
         const Voxel neighbor(get_neighbor(voxel, rnd));
         boost::shared_ptr<const VoxelPool> mt(neighbor.get_voxel_pool());
-        const std::string serial(mt->is_vacant() ? "" : mt->species().serial());
-        if (serial == loc)
+        if (mt->species().serial() == loc)
         {
             tmp.push_back(neighbor);
         }

From 0f3e1c4211b71ce768a45be91d6a12dfe1d59be6 Mon Sep 17 00:00:00 2001
From: Suguru Kato <salmon.bioinformatics@gmail.com>
Date: Thu, 23 Jan 2020 16:07:05 +0900
Subject: [PATCH 21/21] :bug: Fix reactions with OffLatticeSpace

---
 ecell4/spatiocyte/SpatiocyteReactions.cpp |  8 +++----
 ecell4/spatiocyte/SpatiocyteWorld.cpp     | 29 ++++++++++++++++-------
 ecell4/spatiocyte/SpatiocyteWorld.hpp     |  3 +--
 ecell4/spatiocyte/StepEvent.cpp           |  8 ++-----
 4 files changed, 27 insertions(+), 21 deletions(-)

diff --git a/ecell4/spatiocyte/SpatiocyteReactions.cpp b/ecell4/spatiocyte/SpatiocyteReactions.cpp
index 95dc3d7f8..9395aed9d 100644
--- 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
index e286be8f7..af30ceccd 100644
--- 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
index 52c10c779..ea99e2d8d 100644
--- 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
index 65e0c2fd0..585e55b2d 100644
--- 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());