Repository: projectmesa/mesa
Pull request #1994 (2024-01-23 20:51:29)

Add cell-centric discrete spaces (experimental)
## Summary
This PR introduces an alteranative implementation for discrete spaces. This implementation centers on the explicit inclusion of a Cell class. Agents can occupy cells. Cells have connections, specifying their neighbors. The resulting classes offer a cell centric API where agents interact with a cell, and query the cell for its neighbors. 

To capture a collection of cells, and their content (_i.e._, Agents), this PR adds a new CellCollection class. This is an immutable collection of cells with convenient attribute accessors to the cells, or their agents. 

This PR also includes a CellAgent class which extends the default Agent class by adding a `move_to` method that works in conjunction with the new discrete spaces. 

From a performance point of view, the current code is a bit slower in building the grid and cell data structure, but in most use cases this increase in time for model initialization will be more than offset by the faster retrieval of neighboring cells and the agents that occupy them.

## Motive
The PR emerged out of various experiments aimed at improving the performance of the current discrete space code. Moreover, it turned out that a cell centric API resolved various open issues (_e.g._, #1900, #1903, #1953). 


## Implementation
The key idea is to have Cells with connections, and using this to generate neighborhoods for a given radius. So all discrete space classes are in essence a [linked data structure](https://en.wikipedia.org/wiki/Linked_data_structure).

The cell centric API idea is used to implement 4 key discrete space classes: OrthogonalMooreGrid, OrthogonalVonNeumannGrid (alternative for SingleGrid and MultiGrid, and moore and von Neumann neighborhood) , HexGrid (alternative for SingleHexGrid and MultiHexGrid), and Network (alternative for NetworkGrid). Cells have a capacity, so there is no longer a need for seperating Single and Multi grids. Moore and von Neumann reflect different neighborhood connections and so are now implemented as seperate classes. 

## Usage examples
Below are a few examples demonstrating the functionality of the cell spaces.

#### Creating a Grid and Adding Agents
```python
from mesa import Model
from mesa.experimental.cell_space import OrthogonalMooreGrid, CellAgent

class MyAgent(CellAgent):
    def __init__(self, unique_id, model):
        super().__init__(unique_id, model)

# Initialize a model and grid
model = Model()
grid = OrthogonalMooreGrid([10, 10], torus=True, capacity=1, random=model.random)

# Create an agent and add it to a cell
agent = MyAgent(1, model)
cell = grid.select_random_empty_cell()
agent.move_to(cell)
```

#### Agent Movement
```python
# Assuming `agent` and `grid` are already defined
# Move the agent to a random neighboring cell
new_cell = agent.cell.neighborhood().select_random_cell()
agent.move_to(new_cell)
```

#### Querying the Neighborhood
```python
# Get all neighboring agents within a radius of 1
neighbors = agent.cell.neighborhood(radius=1).agents

# Iterate over neighbors and perform actions
for neighbor in neighbors:
    # Example action: print neighbor's ID
    print(neighbor.unique_id)
```

#### Agent Reproduction
```python
# Agent reproduces, creating a new agent in the same cell
def reproduce(agent):
    offspring = MyAgent(model.next_id(), model)
    offspring.move_to(agent.cell)
    model.schedule.add(offspring)

# Call reproduce function for an agent
reproduce(agent)
```

#### Removing an Agent
```python
# Remove an agent from the simulation
def remove_agent(agent):
    agent.cell.remove_agent(agent)  # Remove from cell
    model.schedule.remove(agent)    # Remove from the scheduler

# Call remove_agent function for an agent
remove_agent(agent)
```

#### Checking Cell Properties
```python
# Check if a cell is empty or full before moving an agent
if cell.is_empty:
    agent.move_to(cell)
elif cell.is_full:
    print("Cell is full. Agent cannot move here.")
```

EwoutH: There’s _a lot_ in here that I like!

let me know when you want any input / a review.
EwoutH: If you add `matplotlib` here, the current error in the benchmarks is resolved.

https://github.com/projectmesa/mesa/blob/7b0d3f12f0fb2a0f2dc3ec1bcca4a95a82ff6cee/.github/workflows/benchmarks.yml#L31
Corvince: > If you add `matplotlib` here, the current error in the benchmarks is resolved.
> 
> https://github.com/projectmesa/mesa/blob/7b0d3f12f0fb2a0f2dc3ec1bcca4a95a82ff6cee/.github/workflows/benchmarks.yml#L31

Thanks, I now tried to bypass the jupyter stuff, but I am not sure if this is working. If I change the benchmark.yml won't the action stop working? Because of the pull_request_target thingy?
Corvince: /rerun benchmarks
EwoutH: > If I change the benchmark.yml won't the action stop working?

Yeah it might, honestly I don't know.
rht: Had 3 comments at https://github.com/projectmesa/mesa/commit/2e78940728f2b946abcdca571fd4e85343849a50#comments.
rht: /rerun-benchmarks
github-actions[bot]: Performance benchmarks:

| Model        | Size   | Init time [95% CI]     | Run time [95% CI]      |
|--------------|--------|------------------------|------------------------|
| Schelling    | small  | 🔵 +0.0% [-0.3%, +0.4%] | 🔵 -0.1% [-0.3%, +0.1%] |
| Schelling    | large  | 🔵 -0.5% [-1.1%, +0.1%] | 🔵 -3.9% [-6.4%, -1.2%] |
| WolfSheep    | small  | 🔵 -2.4% [-2.8%, -1.9%] | 🔵 -0.3% [-0.4%, -0.1%] |
| WolfSheep    | large  | 🔵 -2.2% [-3.0%, -1.5%] | 🔵 -3.5% [-4.4%, -2.3%] |
| BoidFlockers | small  | 🔵 -1.7% [-2.2%, -1.1%] | 🔵 +0.4% [-0.0%, +0.9%] |
| BoidFlockers | large  | 🔵 +0.2% [-0.4%, +0.9%] | 🔵 -0.5% [-1.0%, -0.0%] |

Corvince: @EwoutH it looks like the benchmark wasn't run for this PR after the comment trigger. You can see that from the action. Maybe you can look into that?
EwoutH: >  it looks like the benchmark wasn't run for this PR after the comment trigger.

You forgot the dash between rerun and benchmarks. It has to be the exact string (see rht's one).
Corvince: > > it looks like the benchmark wasn't run for this PR after the comment trigger.
> 
> You forgot the dash between rerun and benchmarks. It has to be the exact string (see rht's one).

Ha, I see how my comment was completely miswritten. I already realized after @rht made the correct comment. But in the action that was actually triggered this PR isn't checked out, hence no difference in performance. Performance should actually be about 20% faster.

Here is the exact line of the Action. It should check out this PR, but instead creates a new branch "main"
https://github.com/projectmesa/mesa/actions/runs/7632442972/job/20792667533#step:9:149
quaquel: Thanks for this! I will take a deeper dive at this, hopefully, tonight or otherwise over the weekend. I quickly checked the code this morning and have a few additional comments and suggestions in addition to the various points already raised

I will also move my versions of Schelling and WolfSheep into this (do I have write access)? So we can properly see the change in performance.
EwoutH: > It should check out this PR, but instead creates a new branch "main"

I tried some stuff. Let's see if it works. /rerun-benchmarks
EwoutH: Here we go: /rerun-benchmarks
github-actions[bot]: Performance benchmarks:

| Model        | Size   | Init time [95% CI]      | Run time [95% CI]       |
|--------------|--------|-------------------------|-------------------------|
| Schelling    | small  | 🔵 +1.5% [+0.3%, +2.7%]  | 🔵 +0.8% [+0.1%, +1.5%]  |
| Schelling    | large  | 🔵 -0.4% [-2.0%, +0.9%]  | 🔵 -1.7% [-3.1%, -0.2%]  |
| WolfSheep    | small  | 🔵 +1.6% [+0.5%, +2.6%]  | 🔵 +1.9% [+1.2%, +2.7%]  |
| WolfSheep    | large  | 🔵 +6.4% [+1.7%, +11.7%] | 🔴 +8.2% [+5.3%, +11.1%] |
| BoidFlockers | small  | 🔵 -0.8% [-2.5%, +0.8%]  | 🔵 -0.5% [-1.4%, +0.4%]  |
| BoidFlockers | large  | 🔵 +0.1% [-3.0%, +3.3%]  | 🔵 -0.5% [-1.5%, +0.4%]  |

EwoutH: Doesn't work unfortunately. Will research further at a later point. For now, just do the "mark as ready for review" trick.
Corvince: 
> I will also move my versions of Schelling and WolfSheep into this (do I have write access)? So we can properly see the change in performance.

I invited you as a collaborator just now, so you should have write access. Schelling and WolfSheep are already part of this PR, but the benchmarks are currently not working because I moved the CellSpace into experimental. That means on importing cell_space (for the benchmarks) jupyter_viz is also imported and this leads to a unresolved dependency. Quick solution would be to move cell_space into mesas main directory. Then we can properly retrigger the benchmarks and we should finally see the difference. 


rht: I merged #2001. /rerun-benchmarks
quaquel: I have made some choices in light of all the comments and feedback. Feel free to discuss if you disagree. In short:

* added NetworkGrid
* Renamed various classes
* Added an intermediate Grid class from which OrthogonalGrid and HexGrid derive
* Added a CellAgent class. This isolates this experimental code from the rest of MESA. I also added a `move_to` method to this agent to handle movement.
* I changed the handling of empties closer to the old code. This does introduce coupling between a Cell and its owning "Space" with the various drawbacks that it entails. However, `select_random_empty_cell` now works across OrthogoonalGrid, HexGrid, and NetworkGrid.

compared to my code from last week, this is still a bit slower. I'll try to investigate later. 
quaquel: > compared to my code from last week, this is still a bit slower. I'll try to investigate later.

I fixed the performance issues. The code base is mainly complete. @EwoutH, if possible, let's take 5 minutes tomorrow to discuss property grids and how they might be made to work with this. It's not on the critical path (just as with named connections), but if there is a simple solution, it would be nice to have.

random number generation is still a problem. Once #1981 is addressed, we should use that solution here. 

I still want to add at least some tests before even thinking about merging this. 


Corvince: > I have made some choices in light of all the comments and feedback. Feel free to discuss if you disagree. In short:
> 

Thanks for all the changes!

> * Added a CellAgent class. This isolates this experimental code from the rest of MESA. I also added a `move_to` method to this agent to handle movement.

I don't think this isolates it too much, this is something that has been on several wish lists. But if we manage to find an appropriate implementation for ContinousSpace that is also somewhat "cell-based" (im working on this issue) , these changes could also be propagated to the basic Agent. I know there are ABMs without a space component, but its cell-centric methods don't need to be used. and I don't think there will be too many.

However, I wonder how we will go along with distinguishing cells from positions, esp. in such methods. Currently it only supports cells directly, but you would have to go through model.space to get a reference. I guess thats another point for coupling cells to their owners space. 

> * I changed the handling of empties closer to the old code. This does introduce coupling between a Cell and its owning "Space" with the various drawbacks that it entails. However, `select_random_empty_cell` now works across OrthogoonalGrid, HexGrid, and NetworkGrid.

I still have hopes that we don't need this coupling. I think we can use raise your select_random_empty_cell function from NetworkGrid to the space base class, because it should work for all spaces. Then we can have more optimized implementations for the grids. But re: coupling. I don't like that the cells manipulates the owners attributes (_empties) directly. This is already based on some assumptions how the owner is implemented. This is basically a performance optimization, but those shouldn't dictate how we implement our methods. But yes, without it, it can be quite slow. Not sure how to proceed with this one. 

quaquel: > I don't think this isolates it too much, this is something that has been on several wish lists.

I just meant that it keeps all experimental code within a single location. That is, in this way, we don't touch any legacy code and thus don't have to worry about backward compatibility and breaking existing features. We can simply go to 2.4 with this as new experimental code.

I look forward to seeing your work on continous spaces and I agree that having some space based movement methods within the base agent class is not a problem perse. 

> I still have hopes that we don't need this coupling. I think we can use raise your select_random_empty_cell function from NetworkGrid to the space base class, because it should work for all spaces. 

Will do.

> But re: coupling. I don't like that the cells manipulates the owners attributes (_empties) directly. This is already based on some assumptions how the owner is implemented. This is basically a performance optimization, but those shouldn't dictate how we implement our methods. 

Conceptually, I agree. Within the current structure I don't really see another good solution. Some possible directions I still see are (1) to track the number of tries and use this to estimate when we have to construct the full collection of empty cells; (2) base it on the number of agents in the model, but this introduces tight coupling between model and space; (3) pub/sup where Cells fire events when changing their occupancy. Spaces can then subscribe if needed. 
Corvince:  
> Conceptually, I agree. Within the current structure I don't really see another good solution. Some possible directions I still see are (1) to track the number of tries and use this to estimate when we have to construct the full collection of empty cells; (2) base it on the number of agents in the model, but this introduces tight coupling between model and space; (3) pub/sup where Cells fire events when changing their occupancy. Spaces can then subscribe if needed.

If you agree conceptually I think we shouldn't already break this concept at such an early stage. I tried 1) and think thats a good intermediate solution. Just didn't had the time to find good values. 2) is similar to what we already have, but that actually doesn't always work, because in a MultiGrid you can have lots of agents and still lots of empty cells (this is also problematic in the legacy implementation). 3) sounds like the best solution so far, but requires some work and depends on a general approach towards event sourcing. 

Again, in my opinion we shouldn't break good concepts yet just for the sake of performance optimizations. Until the grid is starting to get full an unoptimized version of 1) already leads to good performance.  

quaquel: > If you agree conceptually I think we shouldn't already break this concept at such an early stage

It happened largely because of me trying to chase down the performance difference between this branch and my code from over a week ago. It turned out that this performance difference was unrelated to this code and I did not change it back. 

> I tried 1) and think thats a good intermediate solution. Just didn't had the time to find good values

I guess we can use the current legacy code to create a number. Will take a look once I have time. The current cutoff implies some expected value, so it should not be too hard. 
quaquel: I have started adding tests for this code. Coverage is not complete yet. I have focused on the initialization of OrthogonalGrid, HexGrid, and NetworkGrid, as well as the neighborhood calculation for various radii and types of grids. This has resulted in a few minor changes in the code. I think this code is getting close to the point that it is ready for review. Still on my to-do list:

- [ ] make Cell and DiscreteSpace again independent and change `Grid.get_random_empty_cell` code as discussed above.
- [ ] add tests for CellCollection, Cell, and CellAgent
- [ ] make neighborhood tests more rigorous. At the moment I test explicitly the size and content of `_connections`, but I only test the size of the `neighborhood` for radius 1-3 on OrthogonalGrid (Moore, von Neumann) and HexGrids. I want to check explicitly the content of the return `neighborhood` as well. 
- [ ] Add some tests for directed graphs and NetworkGrid. 

Further down the line, we can worry about 
1. the integration of property grids. I don't think it is overly complicated, but I don't have the bandwidth for it right now
2. make `_connections` a dict with named connections to enable funky space traversal. Again, see the ideas above

Did I miss anything? If so, please let me know.
Corvince: That's great, thank you for your effort! Tomorrow I want to restructure the file, or better to say divide it into different files, because I think it already has become quite unwieldy. Is that okay for you? I don't want to interfere with your efforts 
quaquel: It's 400 lines, so 25% of the current space module. But indeed, it is getting to a point where it is a bit unwieldy. So I am okay with a reorganization. I am not sure, however, how to cut it up cleanly. Also, if possible can you then reorganize the tests as well so they follow whatever new organization you impose?

I won't work on this tomorrow daytime. I might start working on this again around 19:00 CET. Would that work for you?
Corvince: Yes that should work, otherwise I'll postpone it. Thanks for the reminder to also adjust the test structure 
Corvince: > I am not sure, however, how to cut it up cleanly. 

Simple, I cut it up into each class. Its somewhat less used in the Python world for whatever reason, but I really like to have one file per class. I found it makes switching between code points much easier, because you don't loose your scroll position when jumping between classes. But its still all under the "cell_space" namespace, so the models and tests still use the same import structure. I don't think this is something we need to impose, but I would like to keep it this way for this experimental feature. 

The only exception I made to this is the grid.py file. I think from the content it is quite imminent that its three classes can potentially be merged into a single grid. The only real difference between them is the "directions" list, so I think its rather a question of good naming (we only need a parameter that differentiates between moore, von Neumann and hex).

I kept the tests together for now, I think currently it still makes sense to keep them together. 
quaquel: Thanks for this. I already had a quick look and it makes sense to me (although it gives me a bit of a JAVA feel). 

I have to think about merging the three grid classes. I see your point. Only the `_connect_single_cell` is really different, and thus, it might be handled through a keyword argument. My pushback would be that orthogonal and hex grids are conceptually quite different and distinct notions. Moreover, I am typically fine with a slightly richer set of methods/functions/classes if they represent distinct thinks over solving everything through keyword arguments. 
rht: > Simple, I cut it up into each class. Its somewhat less used in the Python world for whatever reason, but I really like to have one file per class. I found it makes switching between code points much easier, because you don't loose your scroll position when jumping between classes. But its still all under the "cell_space" namespace, so the models and tests still use the same import structure. I don't think this is something we need to impose, but I would like to keep it this way for this experimental feature.

Partially disagree with the way you split the classes here. Either:
- the abstract definitions should be in 1 file, but the specific instances (grid, continuous space, network) be in separate file
- the other way around

While I agree with the general direction of splitting up the classes here, but not all:
1. your reasoning is due to your specific IDE choice/configuration/taste. Several IDE's have gone beyond the boundary of files as organization of code, and have class-based navigation in the left sidebar.
2. as you said, this is a Python project, not Java or MATLAB, and should be less surprising to people who are used to Python project organization.
rht: > the abstract definitions should be in 1 file, but the specific instances (grid, continuous space, network) be in separate file

I meant this:
- cell_space.py
- grid.py
- continuous_space.py
- network.py
Not
- cell_space.py
- implementations.py
quaquel: I suggest we keep the current organization while this is experimental. The code base will likely expand to make this work with property grids and when we add names to the connections. Once the code base is stable, we can worry about the exact organization. 

Note: for the user of MESA it does not matter because, from the outside, it appears as a single namespace. For developers, it means clicking between files instead of through a file with some structured outline. Again, it's not something to overly obsess about. 
rht: The code organization concern can be quickly resolved, once for all, for the code convention throughout the codebase. I find it tiring having to keep track of unresolved concerns until much later.

> For developers, it means clicking between files instead of through a file with some structured outline.

For developers using [VS Code](https://stackoverflow.com/questions/67476757/object-browser-class-browser-for-visual-studio-code), [IntelliJ](https://stackoverflow.com/questions/32511392/in-intellij-is-there-any-package-class-browser-like-eclipses-java-browsing-pers), and [Vim/Neovim](https://www.reddit.com/r/neovim/comments/y9g4is/is_there_a_plugin_similar_to_visual_studios/), one can configure the code browser to be class-based browser instead of file-based.

I know this is bikeshedding, but once again, needs to be resolved once for all. And I'd prefer to follow the standard convention if there is such thing like `black`/`gofmt`/PEP/NEP of this.
quaquel: > I know this is bikeshedding, but once again, needs to be resolved once for all. And I'd prefer to follow the standard convention if there is such thing like black/gofmt/PEP/NEP of this.

Happy to discuss this, but can we do so in a separate place so this stays focussed on the PR? Like you, it would be nice to have some agreements (also on imports) that we can just use as a general code style. 

I don't know of, nor could I find anything so far that gives clear guidance. I can only find very general advice:

1. Up to roughly 400 lines of code in a single module
2. group related classes and functions
3. group modules into packages and sub-packages.
Corvince: 
> I know this is bikeshedding, but once again, needs to be resolved once for all. And I'd prefer to follow the standard convention if there is such thing like `black`/`gofmt`/PEP/NEP of this.

Yes it is bikeshedding and therefore I would kindly ask you to have this discussion elsewhere if you think this is an important matter. I'll reply once here but for future discussions please open a dedicated topic. 

> Partially disagree with the way you split the classes here. Either:
> 
> * the abstract definitions should be in 1 file, but the specific instances (grid, continuous space, network) be in separate file
> * the other way around

I think we are not too far away on this point. The current implementation has all the specific instances in seperate files. As I understand it you would only merge Cell, cellcollection and DiscreteSpace into 1 file. But these are distinct features and therefore shouldn't be in 1 file. Cell, CellCollection and DiscreteSpace follow a hierachical structure and are not tightly coupled. This just leaves cell_agent, but this needs to go near agents anyway once this is out of experimental.

> 
> While I agree with the general direction of splitting up the classes here, but not all:
> 
> 1. your reasoning is due to your specific IDE choice/configuration/taste. Several IDE's have gone beyond the boundary of files as organization of code, and have class-based navigation in the left sidebar.

Its the other way around. All IDEs offer file based structure, some (most) offer class-based navigation. So file-based is more general and less dependent on my specific IDE choice (which I never mentioned).

> 2. as you said, this is a Python project, not Java or MATLAB, and should be less surprising to people who are used to Python project organization.

I think code organisation and language choice are orthogonal to each other, one shouldn't depend on the other. There have been historical reasons for keeping the number of files small, but if you look at modern libraries they tend to spread out their features into many small files and/or group their files per feature. 


quaquel: I have added additional tests. Coverage is pretty complete, but I want to add a few more neighborhood tests just to be sure it's all working correctly. 

I also added two keyword arguments. First, random is now passed around correctly. So, all classes, in principle, use the same random number generator. As long as the DiscreteSpace is initialized with the models rng, everything else works fine. 

I also added `CellKlass` as an optional argument to DiscreteSpace and its subclasses. This is a form of dependency injection. It also makes it trivial for users to subclass Cell while still using MESA predefined DiscreteSpaces.  If, for example, you want Cells with Agent-like behavior, this is now easy to do. 

I still need time to look carefully at the `select_random_empty_cell` code. I have, however, removed already the tight coupling.
quaquel: I have added all the docstrings @EwoutH requested. 
Corvince: I know we are already late in the review process, but I just added two major updates. Both relate to Grids.

1) I added an optional `neighbord_func` argument to the init function. This function can take a cell as its argument and derive the neighborhood directions for it. This makes it a lot more flexible with respect to custom neighborhoods. Based on this change I split OrthogonalGrid into OrthogonalMooreGrid and OrthogonalVonNeumannGrid. This simplifies the code quite a bit, because they now just define their own default neighborhood_func. 

2) I modified Grid to allow n-dimensional grids! API-wise the only, but major change is that instead of `Grid(width, height, ...)` we now have `Grid(dimensions, ...)`. So for the examples I had to change `Grid(width, height, ...)` to `Grid([width, height], ...)`. The internal code changes are quite minimal, but of course this opens up a lot of possibilities. 


This should be up for discussion. I think the code can be improved a bit, maybe even automatically calculate n-dimensional moore/von-Neumann neighborhoods. But I didn't want to go this route further until we agree to keep those changes. 
quaquel: I already had a quick look. 

1. I am not sure I fully understand `neighborhood_func`. Why make it an optional keyword argument if it could be handled through subclassing? I am fine with the Moore vs. von Neuman split into separate classes, which makes good sense to me. It is more explicit this way.
2. I like this direction, and the use of `dimensions` makes sense. Does it already work for, say, 3-D grids? 
Corvince: 1. I actually am not a big fan of oop and for me it's more natural to just define a custom function and pass that to Grid. But that was just the start of my journey to the implementation. But looking at how it is now I fully agree that it should just be handled through subclassing, which feels more natural and less out of place in this context. 

2. I honestly haven't tested it 😅 but it should, the general code is actually quite simple. I just assumed it would work for higher dimensions when it already works for 2 dimensions. The code just iterates over all possible combinations of its dimensions.
quaquel: 1. I broadly agree that composition is often preferable over inheritance, and deep inheritance trees or extensive use of multiple inheritance are a recipe for disaster. But in this case, inheritance is indeed the obvious way to go. 
2. Cell._connections is currently a list. So, one thing worth testing is that cells are not repeated in this list across dimensions. For example, in an orthogonal grid (Moore or Neuman), the cell above the current cell will occur twice if you make the connections naively across three dimensions. 
Corvince: Updated the code to only support inheritance, added parameter validation (mainly to keep hex-grid two-dimensional) and updated moore and von Neumann to work in n dimensions. Added some tests for 3 and 4 dimensions. 

Not sure if Moore and von Neumann is the right nomenclature for n-dimensional grids?
quaquel: > Not sure if Moore and von Neumann is the right nomenclature for n-dimensional grids?

according to [wikipedia](https://en.wikipedia.org/wiki/Von_Neumann_neighborhood) it seems possible to go to higher dimensions. Von Neumann is also applicable to 1D, Moore requires at least 2D.
Corvince: Funny, never thought about the implications for 1-D. I added a test for moore-1d. It still works, although as you shown its technically a von-Neumann grid. But at least it doesn't run into any problems when used for 1d. 

I also updated the the space to be truly generic for different Cell classes. 

Finally I messed up the tests and the commit history ;) Will fix this later, but the goal was actually to make the CI happy. 

quaquel: I had a quick look. The problem with the tests is caused by circular references. In both cell_agent.py and cell_collection.py, you import Cell from cell.py, which in turn imports both cell_agent and cell_collection. 

The reason for importing Cell is that you want to define `T = TypeVar("T", bound=Cell)`. Type hinting is quite new to me, so I am not entirely sure how to fix this. Changing this to "T = TypeVar("T", bound="Cell")" seems to fix the problem. I am however not sure if using a string in this context is allowed. It is in other contexts of type hinting. 

Once done, there is one more test that fails related to Network. I'll try to take a closer look after dinner.
Corvince: Seems to be finally fixed and all checks passing!

The type system of Python is indeed... interesting. I work with TypeScript all day with no problems, but in Python... it sometimes feels really awkward. But you were right, "Cell" instead of Cell did the trick and just required the special TYPE_CHECKING imports
quaquel: I guess this is basically ready to be merged then or not?

Type hinting is just weirdly bolted onto the language over the last few releases. It's getting better and better however. Coming from Java myself, I first liked the dynamic nature of Python, but I more and more long for static typing again. 




Corvince: /rerun-benchmarks
EwoutH: Awesome stuff so far!

I will do final review tomorrow morning. Can I squash when merging?

> /rerun-benchmarks

Unfortunatly, I never got that to work properly, so it was changed to be triggered when a label got added in https://github.com/projectmesa/mesa/pull/2002.
quaquel: > I will do final review tomorrow morning. Can I squash when merging?

There are 2 things left: write a more elaborate commit message when merging and indeed squash and merge. I'll try to write up the message later today. 
github-actions[bot]: Performance benchmarks:

| Model        | Size   | Init time [95% CI]               | Run time [95% CI]          |
|--------------|--------|----------------------------------|----------------------------|
| Schelling    | small  | 🔴 +1768.2% [+1758.4%, +1777.8%] | 🔴 +5.3% [+4.9%, +5.6%]    |
| Schelling    | large  | 🔴 +1308.6% [+1111.5%, +1419.5%] | 🔴 +20.3% [+17.3%, +23.1%] |
| WolfSheep    | small  | 🔴 +303.3% [+301.0%, +305.7%]    | 🟢 -14.7% [-16.0%, -13.6%] |
| WolfSheep    | large  | 🔴 +292.0% [+289.1%, +295.2%]    | 🟢 -20.0% [-21.4%, -18.6%] |
| BoidFlockers | small  | 🔵 +0.5% [-0.2%, +1.2%]          | 🔵 -1.2% [-1.9%, -0.6%]    |
| BoidFlockers | large  | 🔵 +0.6% [-0.0%, +1.2%]          | 🔵 -0.7% [-1.2%, -0.2%]    |

quaquel: draft PR message, feel free to edit

# Summary
This PR introduces an alteranative implementation for discrete spaces. This implementation centers on the explicit inclusion of a Cell class. Agents can occupy cells. Cells have connections, specifying their neighbors. The resulting classes offer a cell centric API where agents interact with a cell, and query the cell for its neighbors. The resulting iplementation is about 20% than the current implementation, has a more expressive API, and is much shorter in terms of lines of code. 

To capture a collection of cells, and their content (_i.e._, Agents), this PR adds a new CellCollection class. This is an immutable collection of cells with convenient attribute accessors to the cells, or their agents. 

This PR also includes a CellAgent class which extends the default Agent class by adding a `move_to` method that works in conjunction with the new discrete spaces. 


# Motive
The PR emerged out of various experiments aimed at improving the performance of the current discrete space code. Moreover, it turned out that a cell centric API resolved various open issues (_e.g._, #1900, #1903, #1953). 


# Implementation
The key idea is to have Cells with connections, and using this to generate neighborhoods for a given radius. So all discrete space classes are in essence a [linked data structure](https://en.wikipedia.org/wiki/Linked_data_structure).

The cell centric API idea is used to implement 4 key discrete space classes: OrthogonalMooreGrid, OrthogonalVonNeumannGrid (alternative for SingleGrid and MultiGrid, and moore and von Neumann neighborhood) , HexGrid (alternative for SingleHexGrid and MultiHexGrid), and Network (alternative for NetworkGrid). Cells have a capacity, so there is no longer a need for seperating Single and Multi grids. Moore and von Neumann reflect different neighborhood connections and so are now implemented as seperate classes.  
EwoutH: Thanks! PR message looks good, I added some usage examples. Could you review if they are correct?

What do we think about the performance? WolfSheep got faster, Schelling slower. I think the 13x to 17x increase in init time for Schelling it quite a lot.

quaquel: I have tried to confirm the results locally and I cannot. Locally, the performance difference for Schelling and WolfSheep is roughly the same. 
github-actions[bot]: Performance benchmarks:

| Model        | Size   | Init time [95% CI]               | Run time [95% CI]          |
|--------------|--------|----------------------------------|----------------------------|
| Schelling    | small  | 🔴 +1789.8% [+1778.9%, +1800.1%] | 🔴 +6.5% [+6.1%, +6.8%]    |
| Schelling    | large  | 🔴 +1317.1% [+1111.6%, +1429.0%] | 🔴 +22.1% [+20.2%, +23.8%] |
| WolfSheep    | small  | 🔴 +302.7% [+301.0%, +304.6%]    | 🟢 -15.4% [-16.6%, -14.2%] |
| WolfSheep    | large  | 🔴 +288.0% [+284.9%, +291.4%]    | 🟢 -20.1% [-21.2%, -18.9%] |
| BoidFlockers | small  | 🔵 +0.1% [-0.6%, +0.9%]          | 🔵 -0.6% [-1.3%, +0.2%]    |
| BoidFlockers | large  | 🔵 +0.8% [+0.1%, +1.5%]          | 🔵 -0.6% [-1.1%, -0.0%]    |

Corvince: For me the local results are pretty much the same as on GitHub
quaquel: Ok, good to know. I'll try to investigate as soon as possible. At a minimum, I want to understand why Schelling is slower in runtime while it was initially about 15-20% faster,  in line with WolfSheep.
quaquel: So, the single biggest loss in performance in the initialization is due to the supporting n-dimensional grids. If a revert back to the older `_connect_single_cell` code, this more than half the initialization time. This is still a factor 2 slower than my old code from a month ago, so there are some other places where we have lost performance.
quaquel: So, changing e.g., `OrthogonalMooreGrid` to the following code below gives me an order of magnitude speedup and we are back to about a factor 4 slower init times compared to main. I have done 2 things: get the offsets only once (hardcoded here for 2D) and _connect_single_cell is likewise 2D only. 

So, the question, at least for init times seems to be: how do we want to handle N-dimensional grids? Writing generic code for N dimensions is an order of magnitude slower than writing dedicated code for 2D. Given that 2D is the most common use case, unless we find a way to speed up the generic code, I guess we have to separate 2D grids from n-dimensional grids. 


```python
class OrthogonalMooreGrid(Grid[T]):
    """Grid where cells are connected to their 8 neighbors.

    Example for two dimensions:
    directions = [
        (-1, -1), (-1, 0), (-1, 1),
        ( 0, -1),          ( 0, 1),
        ( 1, -1), ( 1, 0), ( 1, 1),
    ]
    """

    def __init__(
            self,
            dimensions: Sequence[int],
            torus: bool = False,
            capacity: float | None = None,
            random: Random | None = None,
            cell_klass: type[T] = Cell,

    ) -> None:
        super().__init__(dimensions, torus=torus, capacity=capacity, random=random, cell_klass=cell_klass)
        # offsets = self._calculate_neighborhood_offsets(None)
        # fmt: off
        offsets = [
            (-1, -1), (-1, 0), (-1, 1),
            (0, -1), (0, 1),
            (1, -1), (1, 0), (1, 1),
        ]
        # fmt: on

        for cell in self.all_cells:
            self._connect_single_cell(cell, offsets)

    def _connect_single_cell(self, cell: T, offsets: list[tuple[int, int]]) -> None:
        height, width = self.dimensions

        i, j = cell.coordinate
        for di, dj in offsets:
            ni, nj = (i + di, j + dj)
            if self.torus:
                ni, nj = ni % height, nj % width
            if 0 <= ni < height and 0 <= nj < width:
                cell.connect(self._cells[ni, nj])
```
EwoutH: The current grid also only supports 2D right? For now, and this PR, I would say let's do optimized 2D only. We can always extend the use case further in future.
Corvince: I was worried that performance might not be optimal, but didn't expect (and didn't test) that it would be that bad. 

Still I think it would be a cool feature to keep, but totally agree we should just have an optimized 2D implementation as well
quaquel: I am trying to get to a clean 2d / ND implementation. It will be different across Moore, von Neuman and Hex (which is only defined for 2d anyway. I hope to finish the code in the next hour or so.
github-actions[bot]: Performance benchmarks:

| Model        | Size   | Init time [95% CI]            | Run time [95% CI]          |
|--------------|--------|-------------------------------|----------------------------|
| Schelling    | small  | 🔴 +245.6% [+244.2%, +247.1%] | 🔴 +6.1% [+5.7%, +6.5%]    |
| Schelling    | large  | 🔴 +172.8% [+135.4%, +193.8%] | 🔴 +26.2% [+23.6%, +28.6%] |
| WolfSheep    | small  | 🔴 +71.4% [+69.7%, +73.5%]    | 🟢 -15.1% [-16.2%, -13.9%] |
| WolfSheep    | large  | 🔴 +66.9% [+65.5%, +68.5%]    | 🟢 -16.6% [-18.8%, -14.4%] |
| BoidFlockers | small  | 🔵 +2.7% [+1.9%, +3.4%]       | 🔵 +1.7% [+1.1%, +2.4%]    |
| BoidFlockers | large  | 🔵 +3.5% [+2.8%, +4.1%]       | 🔵 +0.8% [-0.1%, +1.5%]    |

quaquel: The new benchmark figures are quite a bit better. I still don't fully understand why for Schelling the code is slower while in previous tests it was faster than the current implementation. Although one possible explanation is that with the update to AgentSet.shuffle the bottlenecks in the code have shifted quite a bit and in particular percentages have changed a lot as well.
quaquel: So I did a bit more digging. Changing Cell to not use a list if `Cell.capacity` is 1, and updating CellCollection.agent accordingly is an order of magnitude performance difference. 

quaquel: I did some more testing. I can't reproduce my initial performance differences so I guess I made a mistake when I thought there was a massive difference. I also tested a dedicated SingleAgentCell next to the normal cell and updated CellCollection.agents to take advantage of this. Again, this did not produce a massive speedup. So, I suggest we leave the code as is for now and just merge it. I have updated the PR text to better reflect the actual performance figures. 

EwoutH: Agreed, we can always optimize more later if needed.

On thing I'm not sure about is already porting our main benchmarks over. Maybe we could add those models to [mesa-examples](https://github.com/projectmesa/mesa-examples), and leave the benchmarks here as is.
github-actions[bot]: Performance benchmarks:

| Model        | Size   | Init time [95% CI]            | Run time [95% CI]          |
|--------------|--------|-------------------------------|----------------------------|
| Schelling    | small  | 🔴 +240.8% [+239.2%, +242.4%] | 🔴 +7.5% [+7.1%, +7.8%]    |
| Schelling    | large  | 🔴 +170.1% [+132.4%, +190.5%] | 🔴 +26.7% [+24.4%, +28.8%] |
| WolfSheep    | small  | 🔴 +72.9% [+71.1%, +74.5%]    | 🟢 -15.1% [-16.1%, -14.2%] |
| WolfSheep    | large  | 🔴 +61.9% [+44.7%, +71.8%]    | 🟢 -15.8% [-17.9%, -13.6%] |
| BoidFlockers | small  | 🔵 +1.8% [+0.8%, +2.7%]       | 🔵 +0.1% [-0.5%, +0.7%]    |
| BoidFlockers | large  | 🔵 +2.0% [+1.2%, +2.7%]       | 🔵 +0.9% [+0.4%, +1.4%]    |

quaquel: In both my datacollection and devs branch, I have added a simple examples folder within the folder with the new experimental code. I am fine with using either mesa-examples or doing what I have been doing in devs and datacollection. Let me know what you prefer and why.


EwoutH: I would suggest this PR shouldn't change the models in `benchmarks`. I can do those last bits however, I will try to do that and merge it tomorrow.
Corvince: I think an example folder inside the experiments sounds good and like the place to be.


But maybe it would still be nice to have one of the benchmarks use this feature, since I think we all agree that cell space is going to be the recommended space implementation eventually
quaquel: I like @Corvince suggestions. It also allows us to track how the performance of the new grid stuff is affected going forward.
EwoutH: We could have two sets of benchmarks, regular ones and experimental ones
quaquel: I am also fine with adding additional benchmarks, but how would that work out if we add additional experimental stuff? For example, if I add devs, which benchmarks would I edit to showcase their impact?
EwoutH: Thought about it a bit, let's keep it simple. Benchmarks are an internal development tool, not user examples. So let's keep them on the cutting edge.

I'm good with merging (please clean ip or squash). I would like a second aproval before doing so.

Thanks for the thremendous effort!
Corvince: @quaquel can you handle the merging, since I opened the PR? I'm swamped with work until next work and won't be able to do so until then
EwoutH: I can also do it

Patch:
diff --git a/benchmarks/Schelling/schelling.py b/benchmarks/Schelling/schelling.py
index c4b7b01b..c7dd3bf1 100644
--- a/benchmarks/Schelling/schelling.py
+++ b/benchmarks/Schelling/schelling.py
@@ -1,12 +1,14 @@
-import mesa
+from mesa import Model
+from mesa.experimental.cell_space import CellAgent, OrthogonalMooreGrid
+from mesa.time import RandomActivation
 
 
-class SchellingAgent(mesa.Agent):
+class SchellingAgent(CellAgent):
     """
     Schelling segregation agent
     """
 
-    def __init__(self, unique_id, model, agent_type):
+    def __init__(self, unique_id, model, agent_type, radius, homophily):
         """
         Create a new Schelling agent.
         Args:
@@ -16,31 +18,32 @@ class SchellingAgent(mesa.Agent):
         """
         super().__init__(unique_id, model)
         self.type = agent_type
+        self.radius = radius
+        self.homophily = homophily
 
     def step(self):
         similar = 0
-        for neighbor in self.model.grid.iter_neighbors(
-            self.pos, moore=True, radius=self.model.radius
-        ):
+        neighborhood = self.cell.neighborhood(radius=self.radius)
+        for neighbor in neighborhood.agents:
             if neighbor.type == self.type:
                 similar += 1
 
         # If unhappy, move:
-        if similar < self.model.homophily:
-            self.model.grid.move_to_empty(self)
+        if similar < self.homophily:
+            self.move_to(self.model.grid.select_random_empty_cell())
         else:
             self.model.happy += 1
 
 
-class Schelling(mesa.Model):
+class Schelling(Model):
     """
     Model class for the Schelling segregation model.
     """
 
     def __init__(
         self,
-        width=40,
         height=40,
+        width=40,
         homophily=3,
         radius=1,
         density=0.8,
@@ -51,7 +54,7 @@ class Schelling(mesa.Model):
         Create a new Schelling model.
 
         Args:
-            width, height: Size of the space.
+            height, width: Size of the space.
             density: Initial Chance for a cell to populated
             minority_pc: Chances for an agent to be in minority class
             homophily: Minimum number of agents of same class needed to be happy
@@ -59,15 +62,18 @@ class Schelling(mesa.Model):
             seed: Seed for Reproducibility
         """
         super().__init__(seed=seed)
-        self.width = width
         self.height = height
+        self.width = width
         self.density = density
         self.minority_pc = minority_pc
-        self.homophily = homophily
-        self.radius = radius
 
-        self.schedule = mesa.time.RandomActivation(self)
-        self.grid = mesa.space.SingleGrid(width, height, torus=True)
+        self.schedule = RandomActivation(self)
+        self.grid = OrthogonalMooreGrid(
+            [height, width],
+            torus=True,
+            capacity=1,
+            random=self.random,
+        )
 
         self.happy = 0
 
@@ -75,11 +81,13 @@ class Schelling(mesa.Model):
         # We use a grid iterator that returns
         # the coordinates of a cell as well as
         # its contents. (coord_iter)
-        for _, pos in self.grid.coord_iter():
+        for cell in self.grid:
             if self.random.random() < self.density:
                 agent_type = 1 if self.random.random() < self.minority_pc else 0
-                agent = SchellingAgent(self.next_id(), self, agent_type)
-                self.grid.place_agent(agent, pos)
+                agent = SchellingAgent(
+                    self.next_id(), self, agent_type, radius, homophily
+                )
+                agent.move_to(cell)
                 self.schedule.add(agent)
 
     def step(self):
@@ -93,13 +101,12 @@ class Schelling(mesa.Model):
 if __name__ == "__main__":
     import time
 
-    # model = Schelling(seed=15, width=40, height=40, homophily=3, radius=1, density=0.625)
+    # model = Schelling(seed=15, height=40, width=40, homophily=3, radius=1, density=0.625)
     model = Schelling(
-        seed=15, width=100, height=100, homophily=8, radius=2, density=0.8
+        seed=15, height=100, width=100, homophily=8, radius=2, density=0.8
     )
 
     start_time = time.perf_counter()
     for _ in range(100):
         model.step()
-
     print(time.perf_counter() - start_time)
diff --git a/benchmarks/WolfSheep/wolf_sheep.py b/benchmarks/WolfSheep/wolf_sheep.py
index 5d4d9fef..11e07a9e 100644
--- a/benchmarks/WolfSheep/wolf_sheep.py
+++ b/benchmarks/WolfSheep/wolf_sheep.py
@@ -9,41 +9,40 @@ Replication of the model found in NetLogo:
     Northwestern University, Evanston, IL.
 """
 
-import mesa
+import math
 
+from mesa import Model
+from mesa.experimental.cell_space import CellAgent, OrthogonalVonNeumannGrid
+from mesa.time import RandomActivationByType
 
-class Animal(mesa.Agent):
-    def __init__(self, unique_id, model, moore, energy, p_reproduce, energy_from_food):
+
+class Animal(CellAgent):
+    def __init__(self, unique_id, model, energy, p_reproduce, energy_from_food):
         super().__init__(unique_id, model)
         self.energy = energy
         self.p_reproduce = p_reproduce
         self.energy_from_food = energy_from_food
-        self.moore = moore
 
     def random_move(self):
-        next_moves = self.model.grid.get_neighborhood(self.pos, self.moore, True)
-        next_move = self.random.choice(next_moves)
-        # Now move:
-        self.model.grid.move_agent(self, next_move)
+        self.move_to(self.cell.neighborhood().select_random_cell())
 
     def spawn_offspring(self):
         self.energy /= 2
         offspring = self.__class__(
             self.model.next_id(),
             self.model,
-            self.moore,
             self.energy,
             self.p_reproduce,
             self.energy_from_food,
         )
-        self.model.grid.place_agent(offspring, self.pos)
+        offspring.move_to(self.cell)
         self.model.schedule.add(offspring)
 
     def feed(self):
         ...
 
     def die(self):
-        self.model.grid.remove_agent(self)
+        self.cell.remove_agent(self)
         self.remove()
 
     def step(self):
@@ -67,8 +66,9 @@ class Sheep(Animal):
 
     def feed(self):
         # If there is grass available, eat it
-        agents = self.model.grid.get_cell_list_contents(self.pos)
-        grass_patch = next(obj for obj in agents if isinstance(obj, GrassPatch))
+        grass_patch = next(
+            obj for obj in self.cell.agents if isinstance(obj, GrassPatch)
+        )
         if grass_patch.fully_grown:
             self.energy += self.energy_from_food
             grass_patch.fully_grown = False
@@ -80,8 +80,7 @@ class Wolf(Animal):
     """
 
     def feed(self):
-        agents = self.model.grid.get_cell_list_contents(self.pos)
-        sheep = [obj for obj in agents if isinstance(obj, Sheep)]
+        sheep = [obj for obj in self.cell.agents if isinstance(obj, Sheep)]
         if len(sheep) > 0:
             sheep_to_eat = self.random.choice(sheep)
             self.energy += self.energy_from_food
@@ -90,7 +89,7 @@ class Wolf(Animal):
             sheep_to_eat.die()
 
 
-class GrassPatch(mesa.Agent):
+class GrassPatch(CellAgent):
     """
     A patch of grass that grows at a fixed rate and it is eaten by sheep
     """
@@ -100,7 +99,7 @@ class GrassPatch(mesa.Agent):
         Creates a new patch of grass
 
         Args:
-            grown: (boolean) Whether the patch of grass is fully grown or not
+            fully_grown: (boolean) Whether the patch of grass is fully grown or not
             countdown: Time for the patch of grass to be fully grown again
         """
         super().__init__(unique_id, model)
@@ -117,7 +116,7 @@ class GrassPatch(mesa.Agent):
                 self.countdown -= 1
 
 
-class WolfSheep(mesa.Model):
+class WolfSheep(Model):
     """
     Wolf-Sheep Predation Model
 
@@ -126,7 +125,6 @@ class WolfSheep(mesa.Model):
 
     def __init__(
         self,
-        seed,
         height,
         width,
         initial_sheep,
@@ -136,7 +134,7 @@ class WolfSheep(mesa.Model):
         grass_regrowth_time,
         wolf_gain_from_food=13,
         sheep_gain_from_food=5,
-        moore=False,
+        seed=None,
     ):
         """
         Create a new Wolf-Sheep model with the given parameters.
@@ -152,6 +150,7 @@ class WolfSheep(mesa.Model):
                                  once it is eaten
             sheep_gain_from_food: Energy sheep gain from grass, if enabled.
             moore:
+            seed
         """
         super().__init__(seed=seed)
         # Set parameters
@@ -161,9 +160,15 @@ class WolfSheep(mesa.Model):
         self.initial_wolves = initial_wolves
         self.grass_regrowth_time = grass_regrowth_time
 
-        self.schedule = mesa.time.RandomActivationByType(self)
-        self.grid = mesa.space.MultiGrid(self.height, self.width, torus=False)
+        self.schedule = RandomActivationByType(self)
+        self.grid = OrthogonalVonNeumannGrid(
+            [self.height, self.width],
+            torus=False,
+            capacity=math.inf,
+            random=self.random,
+        )
 
+        # Create sheep:
         for _ in range(self.initial_sheep):
             pos = (
                 self.random.randrange(self.width),
@@ -171,14 +176,9 @@ class WolfSheep(mesa.Model):
             )
             energy = self.random.randrange(2 * sheep_gain_from_food)
             sheep = Sheep(
-                self.next_id(),
-                self,
-                moore,
-                energy,
-                sheep_reproduce,
-                sheep_gain_from_food,
+                self.next_id(), self, energy, sheep_reproduce, sheep_gain_from_food
             )
-            self.grid.place_agent(sheep, pos)
+            sheep.move_to(self.grid[pos])
             self.schedule.add(sheep)
 
         # Create wolves
@@ -189,21 +189,21 @@ class WolfSheep(mesa.Model):
             )
             energy = self.random.randrange(2 * wolf_gain_from_food)
             wolf = Wolf(
-                self.next_id(), self, moore, energy, wolf_reproduce, wolf_gain_from_food
+                self.next_id(), self, energy, wolf_reproduce, wolf_gain_from_food
             )
-            self.grid.place_agent(wolf, pos)
+            wolf.move_to(self.grid[pos])
             self.schedule.add(wolf)
 
         # Create grass patches
         possibly_fully_grown = [True, False]
-        for _agent, pos in self.grid.coord_iter():
+        for cell in self.grid:
             fully_grown = self.random.choice(possibly_fully_grown)
             if fully_grown:
                 countdown = self.grass_regrowth_time
             else:
                 countdown = self.random.randrange(self.grass_regrowth_time)
             patch = GrassPatch(self.next_id(), self, fully_grown, countdown)
-            self.grid.place_agent(patch, pos)
+            patch.move_to(cell)
             self.schedule.add(patch)
 
     def step(self):
@@ -213,7 +213,7 @@ class WolfSheep(mesa.Model):
 if __name__ == "__main__":
     import time
 
-    model = WolfSheep(15, 25, 25, 60, 40, 0.2, 0.1, 20)
+    model = WolfSheep(25, 25, 60, 40, 0.2, 0.1, 20, seed=15)
 
     start_time = time.perf_counter()
     for _ in range(100):
diff --git a/mesa/experimental/__init__.py b/mesa/experimental/__init__.py
index c04f2258..961b8762 100644
--- a/mesa/experimental/__init__.py
+++ b/mesa/experimental/__init__.py
@@ -1,1 +1,5 @@
 from .jupyter_viz import JupyterViz, make_text, Slider  # noqa
+from mesa.experimental import cell_space
+
+
+__all__ = ["JupyterViz", "make_text", "Slider", "cell_space"]
diff --git a/mesa/experimental/cell_space/__init__.py b/mesa/experimental/cell_space/__init__.py
new file mode 100644
index 00000000..dce296ae
--- /dev/null
+++ b/mesa/experimental/cell_space/__init__.py
@@ -0,0 +1,23 @@
+from mesa.experimental.cell_space.cell import Cell
+from mesa.experimental.cell_space.cell_agent import CellAgent
+from mesa.experimental.cell_space.cell_collection import CellCollection
+from mesa.experimental.cell_space.discrete_space import DiscreteSpace
+from mesa.experimental.cell_space.grid import (
+    Grid,
+    HexGrid,
+    OrthogonalMooreGrid,
+    OrthogonalVonNeumannGrid,
+)
+from mesa.experimental.cell_space.network import Network
+
+__all__ = [
+    "CellCollection",
+    "Cell",
+    "CellAgent",
+    "DiscreteSpace",
+    "Grid",
+    "HexGrid",
+    "OrthogonalMooreGrid",
+    "OrthogonalVonNeumannGrid",
+    "Network",
+]
diff --git a/mesa/experimental/cell_space/cell.py b/mesa/experimental/cell_space/cell.py
new file mode 100644
index 00000000..55264f68
--- /dev/null
+++ b/mesa/experimental/cell_space/cell.py
@@ -0,0 +1,152 @@
+from __future__ import annotations
+
+from functools import cache
+from random import Random
+from typing import TYPE_CHECKING
+
+from mesa.experimental.cell_space.cell_collection import CellCollection
+
+if TYPE_CHECKING:
+    from mesa.experimental.cell_space.cell_agent import CellAgent
+
+
+class Cell:
+    """The cell represents a position in a discrete space.
+
+    Attributes:
+        coordinate (Tuple[int, int]) : the position of the cell in the discrete space
+        agents (List[Agent]): the agents occupying the cell
+        capacity (int): the maximum number of agents that can simultaneously occupy the cell
+        properties (dict[str, Any]): the properties of the cell
+        random (Random): the random number generator
+
+    """
+
+    __slots__ = [
+        "coordinate",
+        "_connections",
+        "agents",
+        "capacity",
+        "properties",
+        "random",
+    ]
+
+    # def __new__(cls,
+    #     coordinate: tuple[int, ...],
+    #     capacity: float | None = None,
+    #     random: Random | None = None,):
+    #     if capacity != 1:
+    #         return object.__new__(cls)
+    #     else:
+    #         return object.__new__(SingleAgentCell)
+
+    def __init__(
+        self,
+        coordinate: tuple[int, ...],
+        capacity: float | None = None,
+        random: Random | None = None,
+    ) -> None:
+        """ "
+
+        Args:
+            coordinate:
+            capacity (int) : the capacity of the cell. If None, the capacity is infinite
+            random (Random) : the random number generator to use
+
+        """
+        super().__init__()
+        self.coordinate = coordinate
+        self._connections: list[Cell] = []  # TODO: change to CellCollection?
+        self.agents = []  # TODO:: change to AgentSet or weakrefs? (neither is very performant, )
+        self.capacity = capacity
+        self.properties: dict[str, object] = {}
+        self.random = random
+
+    def connect(self, other: Cell) -> None:
+        """Connects this cell to another cell.
+
+        Args:
+            other (Cell): other cell to connect to
+
+        """
+        self._connections.append(other)
+
+    def disconnect(self, other: Cell) -> None:
+        """Disconnects this cell from another cell.
+
+        Args:
+            other (Cell): other cell to remove from connections
+
+        """
+        self._connections.remove(other)
+
+    def add_agent(self, agent: CellAgent) -> None:
+        """Adds an agent to the cell.
+
+        Args:
+            agent (CellAgent): agent to add to this Cell
+
+        """
+        n = len(self.agents)
+
+        if self.capacity and n >= self.capacity:
+            raise Exception(
+                "ERROR: Cell is full"
+            )  # FIXME we need MESA errors or a proper error
+
+        self.agents.append(agent)
+
+    def remove_agent(self, agent: CellAgent) -> None:
+        """Removes an agent from the cell.
+
+        Args:
+            agent (CellAgent): agent to remove from this cell
+
+        """
+        self.agents.remove(agent)
+        agent.cell = None
+
+    @property
+    def is_empty(self) -> bool:
+        """Returns a bool of the contents of a cell."""
+        return len(self.agents) == 0
+
+    @property
+    def is_full(self) -> bool:
+        """Returns a bool of the contents of a cell."""
+        return len(self.agents) == self.capacity
+
+    def __repr__(self):
+        return f"Cell({self.coordinate}, {self.agents})"
+
+    # FIXME: Revisit caching strategy on methods
+    @cache  # noqa: B019
+    def neighborhood(self, radius=1, include_center=False):
+        return CellCollection(
+            self._neighborhood(radius=radius, include_center=include_center),
+            random=self.random,
+        )
+
+    # FIXME: Revisit caching strategy on methods
+    @cache  # noqa: B019
+    def _neighborhood(self, radius=1, include_center=False):
+        # if radius == 0:
+        #     return {self: self.agents}
+        if radius < 1:
+            raise ValueError("radius must be larger than one")
+        if radius == 1:
+            neighborhood = {neighbor: neighbor.agents for neighbor in self._connections}
+            if not include_center:
+                return neighborhood
+            else:
+                neighborhood[self] = self.agents
+                return neighborhood
+        else:
+            neighborhood = {}
+            for neighbor in self._connections:
+                neighborhood.update(
+                    neighbor._neighborhood(radius - 1, include_center=True)
+                )
+            if not include_center:
+                neighborhood.pop(self, None)
+            return neighborhood
diff --git a/mesa/experimental/cell_space/cell_agent.py b/mesa/experimental/cell_space/cell_agent.py
new file mode 100644
index 00000000..abc5155a
--- /dev/null
+++ b/mesa/experimental/cell_space/cell_agent.py
@@ -0,0 +1,37 @@
+from __future__ import annotations
+
+from typing import TYPE_CHECKING
+
+from mesa import Agent, Model
+
+if TYPE_CHECKING:
+    from mesa.experimental.cell_space.cell import Cell
+
+
+class CellAgent(Agent):
+    """Cell Agent is an extension of the Agent class and adds behavior for moving in discrete spaces
+
+
+    Attributes:
+        unique_id (int): A unique identifier for this agent.
+        model (Model): The model instance to which the agent belongs
+        pos: (Position | None): The position of the agent in the space
+        cell: (Cell | None): the cell which the agent occupies
+    """
+
+    def __init__(self, unique_id: int, model: Model) -> None:
+        """
+        Create a new agent.
+
+        Args:
+            unique_id (int): A unique identifier for this agent.
+            model (Model): The model instance in which the agent exists.
+        """
+        super().__init__(unique_id, model)
+        self.cell: Cell | None = None
+
+    def move_to(self, cell) -> None:
+        if self.cell is not None:
+            self.cell.remove_agent(self)
+        self.cell = cell
+        cell.add_agent(self)
diff --git a/mesa/experimental/cell_space/cell_collection.py b/mesa/experimental/cell_space/cell_collection.py
new file mode 100644
index 00000000..9ca36589
--- /dev/null
+++ b/mesa/experimental/cell_space/cell_collection.py
@@ -0,0 +1,81 @@
+from __future__ import annotations
+
+import itertools
+from collections.abc import Iterable, Mapping
+from functools import cached_property
+from random import Random
+from typing import TYPE_CHECKING, Callable, Generic, TypeVar
+
+if TYPE_CHECKING:
+    from mesa.experimental.cell_space.cell import Cell
+    from mesa.experimental.cell_space.cell_agent import CellAgent
+
+T = TypeVar("T", bound="Cell")
+
+
+class CellCollection(Generic[T]):
+    """An immutable collection of cells
+
+    Attributes:
+        cells (List[Cell]): The list of cells this collection represents
+        agents (List[CellAgent]) : List of agents occupying the cells in this collection
+        random (Random) : The random number generator
+
+    """
+
+    def __init__(
+        self,
+        cells: Mapping[T, list[CellAgent]] | Iterable[T],
+        random: Random | None = None,
+    ) -> None:
+        if isinstance(cells, dict):
+            self._cells = cells
+        else:
+            self._cells = {cell: cell.agents for cell in cells}
+
+        #
+        self._capacity: int = next(iter(self._cells.keys())).capacity
+
+        if random is None:
+            random = Random()  # FIXME
+        self.random = random
+
+    def __iter__(self):
+        return iter(self._cells)
+
+    def __getitem__(self, key: T) -> Iterable[CellAgent]:
+        return self._cells[key]
+
+    # @cached_property
+    def __len__(self) -> int:
+        return len(self._cells)
+
+    def __repr__(self):
+        return f"CellCollection({self._cells})"
+
+    @cached_property
+    def cells(self) -> list[T]:
+        return list(self._cells.keys())
+
+    @property
+    def agents(self) -> Iterable[CellAgent]:
+        return itertools.chain.from_iterable(self._cells.values())
+
+    def select_random_cell(self) -> T:
+        return self.random.choice(self.cells)
+
+    def select_random_agent(self) -> CellAgent:
+        return self.random.choice(list(self.agents))
+
+    def select(self, filter_func: Callable[[T], bool] | None = None, n=0):
+        # FIXME: n is not considered
+        if filter_func is None and n == 0:
+            return self
+
+        return CellCollection(
+            {
+                cell: agents
+                for cell, agents in self._cells.items()
+                if filter_func is None or filter_func(cell)
+            }
+        )
diff --git a/mesa/experimental/cell_space/discrete_space.py b/mesa/experimental/cell_space/discrete_space.py
new file mode 100644
index 00000000..d2161c5b
--- /dev/null
+++ b/mesa/experimental/cell_space/discrete_space.py
@@ -0,0 +1,65 @@
+from __future__ import annotations
+
+from functools import cached_property
+from random import Random
+from typing import Generic, TypeVar
+
+from mesa.experimental.cell_space.cell import Cell
+from mesa.experimental.cell_space.cell_collection import CellCollection
+
+T = TypeVar("T", bound=Cell)
+
+
+class DiscreteSpace(Generic[T]):
+    """Base class for all discrete spaces.
+
+    Attributes:
+        capacity (int): The capacity of the cells in the discrete space
+        all_cells (CellCollection): The cells composing the discrete space
+        random (Random): The random number generator
+        cell_klass (Type) : the type of cell class
+        empties (CellCollection) : collecction of all cells that are empty
+
+    """
+
+    def __init__(
+        self,
+        capacity: int | None = None,
+        cell_klass: type[T] = Cell,
+        random: Random | None = None,
+    ):
+        super().__init__()
+        self.capacity = capacity
+        self._cells: dict[tuple[int, ...], T] = {}
+        if random is None:
+            random = Random()  # FIXME should default to default rng from model
+        self.random = random
+        self.cell_klass = cell_klass
+
+        self._empties: dict[tuple[int, ...], None] = {}
+        self._empties_initialized = False
+
+    @property
+    def cutoff_empties(self):
+        return 7.953 * len(self._cells) ** 0.384
+
+    def _connect_single_cell(self, cell: T):
+        ...
+
+    @cached_property
+    def all_cells(self):
+        return CellCollection({cell: cell.agents for cell in self._cells.values()})
+
+    def __iter__(self):
+        return iter(self._cells.values())
+
+    def __getitem__(self, key):
+        return self._cells[key]
+
+    @property
+    def empties(self) -> CellCollection:
+        return self.all_cells.select(lambda cell: cell.is_empty)
+
+    def select_random_empty_cell(self) -> T:
+        """select random empty cell"""
+        return self.random.choice(list(self.empties))
diff --git a/mesa/experimental/cell_space/grid.py b/mesa/experimental/cell_space/grid.py
new file mode 100644
index 00000000..cc4b4b9e
--- /dev/null
+++ b/mesa/experimental/cell_space/grid.py
@@ -0,0 +1,206 @@
+from __future__ import annotations
+
+from collections.abc import Sequence
+from itertools import product
+from random import Random
+from typing import Generic, TypeVar
+
+from mesa.experimental.cell_space import Cell, DiscreteSpace
+
+T = TypeVar("T", bound=Cell)
+
+
+class Grid(DiscreteSpace, Generic[T]):
+    """Base class for all grid classes
+
+    Attributes:
+        dimensions (Sequence[int]): the dimensions of the grid
+        torus (bool): whether the grid is a torus
+        capacity (int): the capacity of a grid cell
+        random (Random): the random number generator
+        _try_random (bool): whether to get empty cell be repeatedly trying random cell
+
+    """
+
+    def __init__(
+        self,
+        dimensions: Sequence[int],
+        torus: bool = False,
+        capacity: float | None = None,
+        random: Random | None = None,
+        cell_klass: type[T] = Cell,
+    ) -> None:
+        super().__init__(capacity=capacity, random=random, cell_klass=cell_klass)
+        self.torus = torus
+        self.dimensions = dimensions
+        self._try_random = True
+        self._ndims = len(dimensions)
+        self._validate_parameters()
+
+        coordinates = product(*(range(dim) for dim in self.dimensions))
+
+        self._cells = {
+            coord: cell_klass(coord, capacity, random=self.random)
+            for coord in coordinates
+        }
+        self._connect_cells()
+
+    def _connect_cells(self) -> None:
+        if self._ndims == 2:
+            self._connect_cells_2d()
+        else:
+            self._connect_cells_nd()
+
+    def _connect_cells_2d(self) -> None:
+        ...
+
+    def _connect_cells_nd(self) -> None:
+        ...
+
+    def _validate_parameters(self):
+        if not all(isinstance(dim, int) and dim > 0 for dim in self.dimensions):
+            raise ValueError("Dimensions must be a list of positive integers.")
+        if not isinstance(self.torus, bool):
+            raise ValueError("Torus must be a boolean.")
+        if self.capacity is not None and not isinstance(self.capacity, (float, int)):
+            raise ValueError("Capacity must be a number or None.")
+
+    def select_random_empty_cell(self) -> T:
+        # FIXME:: currently just a simple boolean to control behavior
+        # FIXME:: basically if grid is close to 99% full, creating empty list can be faster
+        # FIXME:: note however that the old results don't apply because in this implementation
+        # FIXME:: because empties list needs to be rebuild each time
+        # This method is based on Agents.jl's random_empty() implementation. See
+        # https://github.com/JuliaDynamics/Agents.jl/pull/541. For the discussion, see
+        # https://github.com/projectmesa/mesa/issues/1052 and
+        # https://github.com/projectmesa/mesa/pull/1565. The cutoff value provided
+        # is the break-even comparison with the time taken in the else branching point.
+        if self._try_random:
+            while True:
+                cell = self.all_cells.select_random_cell()
+                if cell.is_empty:
+                    return cell
+        else:
+            return super().select_random_empty_cell()
+
+    def _connect_single_cell_nd(self, cell: T, offsets: list[tuple[int, ...]]) -> None:
+        coord = cell.coordinate
+
+        for d_coord in offsets:
+            n_coord = tuple(c + dc for c, dc in zip(coord, d_coord))
+            if self.torus:
+                n_coord = tuple(nc % d for nc, d in zip(n_coord, self.dimensions))
+            if all(0 <= nc < d for nc, d in zip(n_coord, self.dimensions)):
+                cell.connect(self._cells[n_coord])
+
+    def _connect_single_cell_2d(self, cell: T, offsets: list[tuple[int, int]]) -> None:
+        i, j = cell.coordinate
+        height, width = self.dimensions
+
+        for di, dj in offsets:
+            ni, nj = (i + di, j + dj)
+            if self.torus:
+                ni, nj = ni % height, nj % width
+            if 0 <= ni < height and 0 <= nj < width:
+                cell.connect(self._cells[ni, nj])
+
+
+class OrthogonalMooreGrid(Grid[T]):
+    """Grid where cells are connected to their 8 neighbors.
+
+    Example for two dimensions:
+    directions = [
+        (-1, -1), (-1, 0), (-1, 1),
+        ( 0, -1),          ( 0, 1),
+        ( 1, -1), ( 1, 0), ( 1, 1),
+    ]
+    """
+
+    def _connect_cells_2d(self) -> None:
+        # fmt: off
+        offsets = [
+            (-1, -1), (-1, 0), (-1, 1),
+            ( 0, -1),          ( 0, 1),
+            ( 1, -1), ( 1, 0), ( 1, 1),
+        ]
+        # fmt: on
+        height, width = self.dimensions
+
+        for cell in self.all_cells:
+            self._connect_single_cell_2d(cell, offsets)
+
+    def _connect_cells_nd(self) -> None:
+        offsets = list(product([-1, 0, 1], repeat=len(self.dimensions)))
+        offsets.remove((0,) * len(self.dimensions))  # Remove the central cell
+
+        for cell in self.all_cells:
+            self._connect_single_cell_nd(cell, offsets)
+
+
+class OrthogonalVonNeumannGrid(Grid[T]):
+    """Grid where cells are connected to their 4 neighbors.
+
+    Example for two dimensions:
+    directions = [
+                (0, -1),
+        (-1, 0),         ( 1, 0),
+                (0,  1),
+    ]
+    """
+
+    def _connect_cells_2d(self) -> None:
+        # fmt: off
+        offsets = [
+                    (-1, 0),
+            (0, -1),         (0, 1),
+                    ( 1, 0),
+        ]
+        # fmt: on
+        height, width = self.dimensions
+
+        for cell in self.all_cells:
+            self._connect_single_cell_2d(cell, offsets)
+
+    def _connect_cells_nd(self) -> None:
+        offsets = []
+        dimensions = len(self.dimensions)
+        for dim in range(dimensions):
+            for delta in [
+                -1,
+                1,
+            ]:  # Move one step in each direction for the current dimension
+                offset = [0] * dimensions
+                offset[dim] = delta
+                offsets.append(tuple(offset))
+
+        for cell in self.all_cells:
+            self._connect_single_cell_nd(cell, offsets)
+
+
+class HexGrid(Grid[T]):
+    def _connect_cells_2d(self) -> None:
+        # fmt: off
+        even_offsets = [
+                        (-1, -1), (-1, 0),
+                    ( 0, -1),        ( 0, 1),
+                        ( 1, -1), ( 1, 0),
+                ]
+        odd_offsets = [
+                        (-1, 0), (-1, 1),
+                    ( 0, -1),       ( 0, 1),
+                        ( 1, 0), ( 1, 1),
+                ]
+        # fmt: on
+
+        for cell in self.all_cells:
+            i = cell.coordinate[0]
+            offsets = even_offsets if i % 2 == 0 else odd_offsets
+            self._connect_single_cell_2d(cell, offsets=offsets)
+
+    def _connect_cells_nd(self) -> None:
+        raise NotImplementedError("HexGrids are only defined for 2 dimensions")
+
+    def _validate_parameters(self):
+        super()._validate_parameters()
+        if len(self.dimensions) != 2:
+            raise ValueError("HexGrid must have exactly 2 dimensions.")
diff --git a/mesa/experimental/cell_space/network.py b/mesa/experimental/cell_space/network.py
new file mode 100644
index 00000000..57c4d492
--- /dev/null
+++ b/mesa/experimental/cell_space/network.py
@@ -0,0 +1,40 @@
+from random import Random
+from typing import Any, Optional
+
+from mesa.experimental.cell_space.cell import Cell
+from mesa.experimental.cell_space.discrete_space import DiscreteSpace
+
+
+class Network(DiscreteSpace):
+    """A networked discrete space"""
+
+    def __init__(
+        self,
+        G: Any,  # noqa: N803
+        capacity: Optional[int] = None,
+        random: Optional[Random] = None,
+        cell_klass: type[Cell] = Cell,
+    ) -> None:
+        """A Networked grid
+
+        Args:
+            G: a NetworkX Graph instance.
+            capacity (int) : the capacity of the cell
+            random (Random):
+            CellKlass (type[Cell]): The base Cell class to use in the Network
+
+        """
+        super().__init__(capacity=capacity, random=random, cell_klass=cell_klass)
+        self.G = G
+
+        for node_id in self.G.nodes:
+            self._cells[node_id] = self.cell_klass(
+                node_id, capacity, random=self.random
+            )
+
+        for cell in self.all_cells:
+            self._connect_single_cell(cell)
+
+    def _connect_single_cell(self, cell):
+        for node_id in self.G.neighbors(cell.coordinate):
+            cell.connect(self._cells[node_id])


Test patch:
diff --git a/tests/test_cell_space.py b/tests/test_cell_space.py
new file mode 100644
index 00000000..050e39a0
--- /dev/null
+++ b/tests/test_cell_space.py
@@ -0,0 +1,465 @@
+import random
+
+import pytest
+
+from mesa import Model
+from mesa.experimental.cell_space import (
+    Cell,
+    CellAgent,
+    CellCollection,
+    HexGrid,
+    Network,
+    OrthogonalMooreGrid,
+    OrthogonalVonNeumannGrid,
+)
+
+
+def test_orthogonal_grid_neumann():
+    width = 10
+    height = 10
+    grid = OrthogonalVonNeumannGrid((width, height), torus=False, capacity=None)
+
+    assert len(grid._cells) == width * height
+
+    # von neumann neighborhood, torus false, top left corner
+    assert len(grid._cells[(0, 0)]._connections) == 2
+    for connection in grid._cells[(0, 0)]._connections:
+        assert connection.coordinate in {(0, 1), (1, 0)}
+
+    # von neumann neighborhood, torus false, top right corner
+    for connection in grid._cells[(0, width - 1)]._connections:
+        assert connection.coordinate in {(0, width - 2), (1, width - 1)}
+
+    # von neumann neighborhood, torus false, bottom left corner
+    for connection in grid._cells[(height - 1, 0)]._connections:
+        assert connection.coordinate in {(height - 1, 1), (height - 2, 0)}
+
+    # von neumann neighborhood, torus false, bottom right corner
+    for connection in grid._cells[(height - 1, width - 1)]._connections:
+        assert connection.coordinate in {
+            (height - 1, width - 2),
+            (height - 2, width - 1),
+        }
+
+    # von neumann neighborhood middle of grid
+    assert len(grid._cells[(5, 5)]._connections) == 4
+    for connection in grid._cells[(5, 5)]._connections:
+        assert connection.coordinate in {(4, 5), (5, 4), (5, 6), (6, 5)}
+
+    # von neumann neighborhood, torus True, top corner
+    grid = OrthogonalVonNeumannGrid((width, height), torus=True, capacity=None)
+    assert len(grid._cells[(0, 0)]._connections) == 4
+    for connection in grid._cells[(0, 0)]._connections:
+        assert connection.coordinate in {(0, 1), (1, 0), (0, 9), (9, 0)}
+
+    # von neumann neighborhood, torus True, top right corner
+    for connection in grid._cells[(0, width - 1)]._connections:
+        assert connection.coordinate in {(0, 8), (0, 0), (1, 9), (9, 9)}
+
+    # von neumann neighborhood, torus True, bottom left corner
+    for connection in grid._cells[(9, 0)]._connections:
+        assert connection.coordinate in {(9, 1), (9, 9), (0, 0), (8, 0)}
+
+    # von neumann neighborhood, torus True, bottom right corner
+    for connection in grid._cells[(9, 9)]._connections:
+        assert connection.coordinate in {(9, 0), (9, 8), (8, 9), (0, 9)}
+
+
+def test_orthogonal_grid_neumann_3d():
+    width = 10
+    height = 10
+    depth = 10
+    grid = OrthogonalVonNeumannGrid((width, height, depth), torus=False, capacity=None)
+
+    assert len(grid._cells) == width * height * depth
+
+    # von neumann neighborhood, torus false, top left corner
+    assert len(grid._cells[(0, 0, 0)]._connections) == 3
+    for connection in grid._cells[(0, 0, 0)]._connections:
+        assert connection.coordinate in {(0, 0, 1), (0, 1, 0), (1, 0, 0)}
+
+    # von neumann neighborhood, torus false, top right corner
+    for connection in grid._cells[(0, width - 1, 0)]._connections:
+        assert connection.coordinate in {
+            (0, width - 1, 1),
+            (0, width - 2, 0),
+            (1, width - 1, 0),
+        }
+
+    # von neumann neighborhood, torus false, bottom left corner
+    for connection in grid._cells[(height - 1, 0, 0)]._connections:
+        assert connection.coordinate in {
+            (height - 1, 0, 1),
+            (height - 1, 1, 0),
+            (height - 2, 0, 0),
+        }
+
+    # von neumann neighborhood, torus false, bottom right corner
+    for connection in grid._cells[(height - 1, width - 1, 0)]._connections:
+        assert connection.coordinate in {
+            (height - 1, width - 1, 1),
+            (height - 1, width - 2, 0),
+            (height - 2, width - 1, 0),
+        }
+
+    # von neumann neighborhood middle of grid
+    assert len(grid._cells[(5, 5, 5)]._connections) == 6
+    for connection in grid._cells[(5, 5, 5)]._connections:
+        assert connection.coordinate in {
+            (4, 5, 5),
+            (5, 4, 5),
+            (5, 5, 4),
+            (5, 5, 6),
+            (5, 6, 5),
+            (6, 5, 5),
+        }
+
+    # von neumann neighborhood, torus True, top corner
+    grid = OrthogonalVonNeumannGrid((width, height, depth), torus=True, capacity=None)
+    assert len(grid._cells[(0, 0, 0)]._connections) == 6
+    for connection in grid._cells[(0, 0, 0)]._connections:
+        assert connection.coordinate in {
+            (0, 0, 1),
+            (0, 1, 0),
+            (1, 0, 0),
+            (0, 0, 9),
+            (0, 9, 0),
+            (9, 0, 0),
+        }
+
+
+def test_orthogonal_grid_moore():
+    width = 10
+    height = 10
+
+    # Moore neighborhood, torus false, top corner
+    grid = OrthogonalMooreGrid((width, height), torus=False, capacity=None)
+    assert len(grid._cells[(0, 0)]._connections) == 3
+    for connection in grid._cells[(0, 0)]._connections:
+        assert connection.coordinate in {(0, 1), (1, 0), (1, 1)}
+
+    # Moore neighborhood middle of grid
+    assert len(grid._cells[(5, 5)]._connections) == 8
+    for connection in grid._cells[(5, 5)]._connections:
+        # fmt: off
+        assert connection.coordinate in {(4, 4), (4, 5), (4, 6),
+                                         (5, 4),         (5, 6),
+                                         (6, 4), (6, 5), (6, 6)}
+        # fmt: on
+
+    # Moore neighborhood, torus True, top corner
+    grid = OrthogonalMooreGrid([10, 10], torus=True, capacity=None)
+    assert len(grid._cells[(0, 0)]._connections) == 8
+    for connection in grid._cells[(0, 0)]._connections:
+        # fmt: off
+        assert connection.coordinate in {(9, 9), (9, 0), (9, 1),
+                                         (0, 9),         (0, 1),
+                                         (1, 9), (1, 0), (1, 1)}
+        # fmt: on
+
+
+def test_orthogonal_grid_moore_3d():
+    width = 10
+    height = 10
+    depth = 10
+
+    # Moore neighborhood, torus false, top corner
+    grid = OrthogonalMooreGrid((width, height, depth), torus=False, capacity=None)
+    assert len(grid._cells[(0, 0, 0)]._connections) == 7
+    for connection in grid._cells[(0, 0, 0)]._connections:
+        assert connection.coordinate in {
+            (0, 0, 1),
+            (0, 1, 0),
+            (0, 1, 1),
+            (1, 0, 0),
+            (1, 0, 1),
+            (1, 1, 0),
+            (1, 1, 1),
+        }
+
+    # Moore neighborhood middle of grid
+    assert len(grid._cells[(5, 5, 5)]._connections) == 26
+    for connection in grid._cells[(5, 5, 5)]._connections:
+        # fmt: off
+        assert connection.coordinate in {(4, 4, 4), (4, 4, 5), (4, 4, 6), (4, 5, 4), (4, 5, 5), (4, 5, 6), (4, 6, 4), (4, 6, 5), (4, 6, 6),
+                                         (5, 4, 4), (5, 4, 5), (5, 4, 6), (5, 5, 4),             (5, 5, 6), (5, 6, 4), (5, 6, 5), (5, 6, 6),
+                                         (6, 4, 4), (6, 4, 5), (6, 4, 6), (6, 5, 4), (6, 5, 5), (6, 5, 6), (6, 6, 4), (6, 6, 5), (6, 6, 6)}
+        # fmt: on
+
+    # Moore neighborhood, torus True, top corner
+    grid = OrthogonalMooreGrid((width, height, depth), torus=True, capacity=None)
+    assert len(grid._cells[(0, 0, 0)]._connections) == 26
+    for connection in grid._cells[(0, 0, 0)]._connections:
+        # fmt: off
+        assert connection.coordinate in {(9, 9, 9), (9, 9, 0), (9, 9, 1), (9, 0, 9), (9, 0, 0), (9, 0, 1), (9, 1, 9), (9, 1, 0), (9, 1, 1),
+                                         (0, 9, 9), (0, 9, 0), (0, 9, 1), (0, 0, 9),             (0, 0, 1), (0, 1, 9), (0, 1, 0), (0, 1, 1),
+                                         (1, 9, 9), (1, 9, 0), (1, 9, 1), (1, 0, 9), (1, 0, 0), (1, 0, 1), (1, 1, 9), (1, 1, 0), (1, 1, 1)}
+        # fmt: on
+
+
+def test_orthogonal_grid_moore_4d():
+    width = 10
+    height = 10
+    depth = 10
+    time = 10
+
+    # Moore neighborhood, torus false, top corner
+    grid = OrthogonalMooreGrid((width, height, depth, time), torus=False, capacity=None)
+    assert len(grid._cells[(0, 0, 0, 0)]._connections) == 15
+    for connection in grid._cells[(0, 0, 0, 0)]._connections:
+        assert connection.coordinate in {
+            (0, 0, 0, 1),
+            (0, 0, 1, 0),
+            (0, 0, 1, 1),
+            (0, 1, 0, 0),
+            (0, 1, 0, 1),
+            (0, 1, 1, 0),
+            (0, 1, 1, 1),
+            (1, 0, 0, 0),
+            (1, 0, 0, 1),
+            (1, 0, 1, 0),
+            (1, 0, 1, 1),
+            (1, 1, 0, 0),
+            (1, 1, 0, 1),
+            (1, 1, 1, 0),
+            (1, 1, 1, 1),
+        }
+
+    # Moore neighborhood middle of grid
+    assert len(grid._cells[(5, 5, 5, 5)]._connections) == 80
+    for connection in grid._cells[(5, 5, 5, 5)]._connections:
+        # fmt: off
+        assert connection.coordinate in {(4, 4, 4, 4), (4, 4, 4, 5), (4, 4, 4, 6), (4, 4, 5, 4), (4, 4, 5, 5), (4, 4, 5, 6), (4, 4, 6, 4), (4, 4, 6, 5), (4, 4, 6, 6),
+                                         (4, 5, 4, 4), (4, 5, 4, 5), (4, 5, 4, 6), (4, 5, 5, 4), (4, 5, 5, 5), (4, 5, 5, 6), (4, 5, 6, 4), (4, 5, 6, 5), (4, 5, 6, 6),
+                                            (4, 6, 4, 4), (4, 6, 4, 5), (4, 6, 4, 6), (4, 6, 5, 4), (4, 6, 5, 5), (4, 6, 5, 6), (4, 6, 6, 4), (4, 6, 6, 5), (4, 6, 6, 6),
+                                            (5, 4, 4, 4), (5, 4, 4, 5), (5, 4, 4, 6), (5, 4, 5, 4), (5, 4, 5, 5), (5, 4, 5, 6), (5, 4, 6, 4), (5, 4, 6, 5), (5, 4, 6, 6),
+                                            (5, 5, 4, 4), (5, 5, 4, 5), (5, 5, 4, 6), (5, 5, 5, 4),             (5, 5, 5, 6), (5, 5, 6, 4), (5, 5, 6, 5), (5, 5, 6, 6),
+                                            (5, 6, 4, 4), (5, 6, 4, 5), (5, 6, 4, 6), (5, 6, 5, 4), (5, 6, 5, 5), (5, 6, 5, 6), (5, 6, 6, 4), (5, 6, 6, 5), (5, 6, 6, 6),
+                                            (6, 4, 4, 4), (6, 4, 4, 5), (6, 4, 4, 6), (6, 4, 5, 4), (6, 4, 5, 5), (6, 4, 5, 6), (6, 4, 6, 4), (6, 4, 6, 5), (6, 4, 6, 6),
+                                            (6, 5, 4, 4), (6, 5, 4, 5), (6, 5, 4, 6), (6, 5, 5, 4), (6, 5, 5, 5), (6, 5, 5, 6), (6, 5, 6, 4), (6, 5, 6, 5), (6, 5, 6, 6),
+                                            (6, 6, 4, 4), (6, 6, 4, 5), (6, 6, 4, 6), (6, 6, 5, 4), (6, 6, 5, 5), (6, 6, 5, 6), (6, 6, 6, 4), (6, 6, 6, 5), (6, 6, 6, 6)}
+        # fmt: on
+
+
+def test_orthogonal_grid_moore_1d():
+    width = 10
+
+    # Moore neighborhood, torus false, left edge
+    grid = OrthogonalMooreGrid((width,), torus=False, capacity=None)
+    assert len(grid._cells[(0,)]._connections) == 1
+    for connection in grid._cells[(0,)]._connections:
+        assert connection.coordinate in {(1,)}
+
+    # Moore neighborhood middle of grid
+    assert len(grid._cells[(5,)]._connections) == 2
+    for connection in grid._cells[(5,)]._connections:
+        assert connection.coordinate in {(4,), (6,)}
+
+    # Moore neighborhood, torus True, left edge
+    grid = OrthogonalMooreGrid((width,), torus=True, capacity=None)
+    assert len(grid._cells[(0,)]._connections) == 2
+    for connection in grid._cells[(0,)]._connections:
+        assert connection.coordinate in {(1,), (9,)}
+
+
+def test_cell_neighborhood():
+    # orthogonal grid
+
+    ## von Neumann
+    width = 10
+    height = 10
+    grid = OrthogonalVonNeumannGrid((width, height), torus=False, capacity=None)
+    for radius, n in zip(range(1, 4), [2, 5, 9]):
+        neighborhood = grid._cells[(0, 0)].neighborhood(radius=radius)
+        assert len(neighborhood) == n
+
+    ## Moore
+    width = 10
+    height = 10
+    grid = OrthogonalMooreGrid((width, height), torus=False, capacity=None)
+    for radius, n in zip(range(1, 4), [3, 8, 15]):
+        neighborhood = grid._cells[(0, 0)].neighborhood(radius=radius)
+        assert len(neighborhood) == n
+
+    with pytest.raises(ValueError):
+        grid._cells[(0, 0)].neighborhood(radius=0)
+
+    # hexgrid
+    width = 10
+    height = 10
+    grid = HexGrid((width, height), torus=False, capacity=None)
+    for radius, n in zip(range(1, 4), [2, 6, 11]):
+        neighborhood = grid._cells[(0, 0)].neighborhood(radius=radius)
+        assert len(neighborhood) == n
+
+    width = 10
+    height = 10
+    grid = HexGrid((width, height), torus=False, capacity=None)
+    for radius, n in zip(range(1, 4), [5, 10, 17]):
+        neighborhood = grid._cells[(1, 0)].neighborhood(radius=radius)
+        assert len(neighborhood) == n
+
+    # networkgrid
+
+
+def test_hexgrid():
+    width = 10
+    height = 10
+
+    grid = HexGrid((width, height), torus=False)
+    assert len(grid._cells) == width * height
+
+    # first row
+    assert len(grid._cells[(0, 0)]._connections) == 2
+    for connection in grid._cells[(0, 0)]._connections:
+        assert connection.coordinate in {(0, 1), (1, 0)}
+
+    # second row
+    assert len(grid._cells[(1, 0)]._connections) == 5
+    for connection in grid._cells[(1, 0)]._connections:
+        # fmt: off
+        assert connection.coordinate in {(0, 0), (0, 1),
+                                                    (1, 1),
+                                         (2, 0), (2, 1)}
+
+    # middle odd row
+    assert len(grid._cells[(5, 5)]._connections) == 6
+    for connection in grid._cells[(5, 5)]._connections:
+        # fmt: off
+        assert connection.coordinate in {(4, 5), (4, 6),
+                                      (5, 4),       (5, 6),
+                                         (6, 5), (6, 6)}
+
+        # fmt: on
+
+    # middle even row
+    assert len(grid._cells[(4, 4)]._connections) == 6
+    for connection in grid._cells[(4, 4)]._connections:
+        # fmt: off
+        assert connection.coordinate in {(3, 3), (3, 4),
+                                      (4, 3),       (4, 5),
+                                         (5, 3), (5, 4)}
+
+        # fmt: on
+
+    grid = HexGrid((width, height), torus=True)
+    assert len(grid._cells) == width * height
+
+    # first row
+    assert len(grid._cells[(0, 0)]._connections) == 6
+    for connection in grid._cells[(0, 0)]._connections:
+        # fmt: off
+        assert connection.coordinate in {(9, 9), (9, 0),
+                                      (0, 9),       (0, 1),
+                                         (1, 9), (1, 0)}
+
+        # fmt: on
+
+
+def test_networkgrid():
+    import networkx as nx
+
+    n = 10
+    m = 20
+    seed = 42
+    G = nx.gnm_random_graph(n, m, seed=seed)  # noqa: N806
+    grid = Network(G)
+
+    assert len(grid._cells) == n
+
+    for i, cell in grid._cells.items():
+        for connection in cell._connections:
+            assert connection.coordinate in G.neighbors(i)
+
+
+def test_empties_space():
+    import networkx as nx
+
+    n = 10
+    m = 20
+    seed = 42
+    G = nx.gnm_random_graph(n, m, seed=seed)  # noqa: N806
+    grid = Network(G)
+
+    assert len(grid.empties) == n
+
+    model = Model()
+    for i in range(8):
+        grid._cells[i].add_agent(CellAgent(i, model))
+
+    cell = grid.select_random_empty_cell()
+    assert cell.coordinate in {8, 9}
+
+
+def test_cell():
+    cell1 = Cell((1,), capacity=None, random=random.Random())
+    cell2 = Cell((2,), capacity=None, random=random.Random())
+
+    # connect
+    cell1.connect(cell2)
+    assert cell2 in cell1._connections
+
+    # disconnect
+    cell1.disconnect(cell2)
+    assert cell2 not in cell1._connections
+
+    # remove cell not in connections
+    with pytest.raises(ValueError):
+        cell1.disconnect(cell2)
+
+    # add_agent
+    model = Model()
+    agent = CellAgent(1, model)
+
+    cell1.add_agent(agent)
+    assert agent in cell1.agents
+
+    # remove_agent
+    cell1.remove_agent(agent)
+    assert agent not in cell1.agents
+
+    with pytest.raises(ValueError):
+        cell1.remove_agent(agent)
+
+    cell1 = Cell((1,), capacity=1, random=random.Random())
+    cell1.add_agent(CellAgent(1, model))
+    assert cell1.is_full
+
+    with pytest.raises(Exception):
+        cell1.add_agent(CellAgent(2, model))
+
+
+def test_cell_collection():
+    cell1 = Cell((1,), capacity=None, random=random.Random())
+
+    collection = CellCollection({cell1: cell1.agents}, random=random.Random())
+    assert len(collection) == 1
+    assert cell1 in collection
+
+    rng = random.Random()
+    n = 10
+    collection = CellCollection([Cell((i,), random=rng) for i in range(n)], random=rng)
+    assert len(collection) == n
+
+    cell = collection.select_random_cell()
+    assert cell in collection
+
+    cells = collection.cells
+    assert len(cells) == n
+
+    agents = collection.agents
+    assert len(list(agents)) == 0
+
+    cells = collection.cells
+    model = Model()
+    cells[0].add_agent(CellAgent(1, model))
+    cells[3].add_agent(CellAgent(2, model))
+    cells[7].add_agent(CellAgent(3, model))
+    agents = collection.agents
+    assert len(list(agents)) == 3
+
+    agent = collection.select_random_agent()
+    assert agent in set(collection.agents)
+
+    agents = collection[cells[0]]
+    assert agents == cells[0].agents