Repository: projectmesa/mesa
Pull request #1898 (2023-12-04 12:36:43)

space: Implement PropertyLayer and _PropertyGrid
### Summary
This PR introduces a system for managing cell properties within Mesa's grid classes. It adds a `PropertyLayer` class and extends the `SingleGrid` and `MultiGrid` classes to support cell properties, enhancing the capability of Mesa for more detailed environmental modeling in simulations.

### Motive
The addition of cell properties brings a new dimension to agent-based models, allowing them to simulate environmental factors (e.g., terrain types, resources, pollution levels) more realistically. 

From a technical perspective, cell properties reduce computational overhead by negating the need for numerous agent objects for static environmental features. Conceptually, they provide a clearer and more intuitive representation of the environment, streamlining model complexity and facilitating analysis.

This feature enhances Mesa's modeling capabilities by enabling dynamic interaction between agents and their environment. It allows for simulations where environmental conditions directly influence agent behavior, thus broadening the scope and depth of possible models.

### Implementation
![pako_eNrdVd9vmzAQ_lcsP9GNROR3gqpK0yrtZZUm5W1CQi6-EKtgI2Pasiz_-w5CgqGk27S38ZDA3Xd3n-8-2wcaKQ7Up1HC8vxesFizNJAEn9pCvmmVgTblV1aCJoeTq3pGW6OFjIlkKVhWIQ15Edzse7Y9iHhvLKMs0qwcS860ZiXhzLDW-TEHE0aQJM5npTQXkhkgmcqFEUq65JklB](https://github.com/projectmesa/mesa/assets/15776622/456d2fca-699c-4c9d-82cc-dd0a2965465f)

#### PropertyLayer Class
The `PropertyLayer` class represents a layer of properties associated with the cells of a grid. It includes the following functionalities:

- **Initialization:** The constructor initializes a 2D NumPy array to store property values for each cell, with a specified default value and data type.
- **Cell Value Management:** Methods to set or modify the value of individual cells or batches of cells, either unconditionally or based on specified conditions.
- **Property Manipulation:** Functions to apply operations (such as lambda functions or NumPy universal functions) to modify cell values, either individually or in a batch mode.
- **Property Querying:** Capabilities to select cells based on certain conditions and perform aggregate operations across all cells (e.g., sum, average).

#### _PropertyGrid Class
The `_PropertyGrid` class extends the existing `_Grid` class, adding support for managing multiple `PropertyLayer` instances. Key features include:

- **Property Layer Management:** Methods to add and remove property layers from the grid. Each property layer is identified by a unique name and is managed independently.
- **Interaction with Agents:** Integration of property layers with agent actions, allowing agents to interact with and modify cell properties as part of their behavior.
- **Enhanced Grid Functionalities:** Extension of existing grid functionalities to incorporate property-based conditions. For example, enabling agents to move to cells with specific properties or query neighboring cells based on property values.

#### Integration with Grid Classes
The existing grid classes, such as `SingleGrid` and `MultiGrid`, are derived from `_PropertyGrid`. This inheritance allows them to leverage the property layer functionalities seamlessly. These enhanced grid classes enable modelers to define and manipulate cell properties directly, enriching the modeling capabilities of Mesa.

#### Implementation Details
- The implementation is compatible with existing Mesa models, ensuring no disruption to current users.
- It supports lambda functions and NumPy ufuncs, providing flexibility in defining cell operations.

### API examples
Here are some code examples of how these features could be utilized.

Initializing a Property Layer and Setting Cell Values
```python
import numpy as np
from mesa.space import SingleGrid, PropertyLayer

# Initialize a SingleGrid with a PropertyLayer
grid = SingleGrid(10, 10, False)
property_layer = PropertyLayer("elevation", 10, 10, default_value=0)
grid.add_property_layer(property_layer)

# Set a specific cell's value
grid.properties["elevation"].set_cell((5, 5), 100)

# Set multiple cells' values based on a condition
condition = lambda x: x < 50  # Condition to check each cell's value
grid.properties["elevation"].set_cells(60, condition)
```
Modifying Cell Values Using Operations.
Note: If performance is important, NumPy [Universal functions](https://numpy.org/doc/stable/reference/ufuncs.html) (ufunc) are often faster.
```python
# Modify a single cell's value with a lambda function
grid.properties["elevation"].modify_cell((3, 3), lambda x: x + 20)

# Modify multiple cells' values with a NumPy ufunc
grid.properties["elevation"].modify_cells(np.multiply, 2)

```
Conditional Operations Based on Cell Properties
```python
# Define a condition function
high_pollution = lambda x: x > 3.0

# Select cells that meet the condition
high_pollution_cells = grid.properties['pollution'].select_cells(condition=high_pollution)

# Reduce pollution in highly polluted cells using a lambda function
reduce_pollution = lambda x: x - 2.0 if x > 3.0 else x
grid.properties['pollution'].modify_cells(operation=reduce_pollution)
```
Aggregating Property Values
We can perform an aggregate operation, such as calculating the sum or mean of all values in the property layer.
```python
# Calculate the total pollution level
total_pollution = grid.properties['pollution'].aggregate_property(np.sum)

# Calculate the average pollution level
average_pollution = grid.properties['pollution'].aggregate_property(np.mean)
```
Complex Cell Selection Based on Condition
We can select cells that meet a complex condition, such as finding cells within a specific range of values.
```python
# Select cells where the value is between 5 and 15
selected_cells = layer.select_cells(lambda x: (x > 5) & (x < 15))

# Print the coordinates of the selected cells
print("Selected Cells:", selected_cells)
```

### Model examples
##### Moving an Agent to a Cell with the Highest Value
In this example, we'll move an agent to the cell with the highest value in a specified property layer.
```Python
from mesa import Agent, Model
from mesa.space import SingleGrid, PropertyLayer
import numpy as np

class MyAgent(Agent):
    def move_to_highest(self):
        # Move the agent to the cell with the highest value in 'temperature' layer
        target_cells = self.model.grid.select_cells(extreme_values={"temperature": "highest"})
        self.model.grid.move_agent_to_one_of(self, target_cells)

class MyModel(Model):
    def __init__(self):
        self.grid = SingleGrid(10, 10, False)
        temperature_layer = PropertyLayer("temperature", 10, 10, default_value=20, dtype=int)
        self.grid.add_property_layer(temperature_layer)

        # Randomly increase temperature in some cells
        for _ in range(30):
            x, y = self.random.randint(0, 9), self.random.randint(0, 9)
            temperature_layer.set_cell((x, y), self.random.randint(20, 40))

        # Add an agent
        self.agent = MyAgent(0, self)
        self.grid.place_agent(self.agent, (5, 5))

# Create model and move agent
model = MyModel()
model.agent.move_to_highest()

print(f"Agent moved to: {model.agent.pos}")
```
##### Example 2: Moving an Agent Based on Multiple Property Conditions
Here, we move an agent to a random cell that satisfies multiple conditions based on different property layers.

```Python
from mesa import Agent, Model
from mesa.space import SingleGrid, PropertyLayer
import numpy as np

class MyAgent(Agent):
    def move_to_optimal_cell(self):
        conditions = {
            "temperature": lambda x: x < 30,
            "humidity": lambda x: x > 50
        }
        target_cells = self.model.grid.select_cells(conditions)
        if target_cells:
            self.model.grid.move_agent_to_one_of(self, target_cells)

class MyModel(Model):
    def __init__(self):
        self.grid = SingleGrid(10, 10, False)
        self.grid.add_property_layer(PropertyLayer("temperature", 10, 10, 20, dtype=int))
        self.grid.add_property_layer(PropertyLayer("humidity", 10, 10, 60.0, dtype=float))

        # Add an agent
        self.agent = MyAgent(0, self)
        self.grid.place_agent(self.agent, (5, 5))

# Create model and move agent
model = MyModel()
model.agent.move_to_optimal_cell()

print(f"Agent moved to: {model.agent.pos}")
```
### Scenario: Forest Fire Simulation

In this simulation, we have a forest with varying levels of temperature, humidity, and elevation. Agents represent firefighters trying to contain a forest fire. The goal is for the firefighters to move towards areas with the highest temperature (indicative of fire) but also consider safety constraints such as avoiding extremely low humidity areas (high fire risk), inaccessible high elevation zones or extremely high temperature zones. They also have a maximum movement radius of 3 (using the neighbor mask) and won't move to a cell where already a firefighter is (using the empty mask).

```python
import numpy as np

from mesa import Agent, Model
from mesa.space import PropertyLayer, SingleGrid
from mesa.time import RandomActivation


class FirefighterAgent(Agent):
    def step(self):
        current_location = self.pos

        # Move to the hottest eligible cell
        self.move_to_hottest_eligible_cell()

        current_temperature = self.model.grid.properties["temperature"].data[self.pos]

        # Decrease the temperature (in Kelvin) by 20% at the agent's current cell
        self.reduce_temperature()

        print(
            f"Firefighter {self.unique_id} moved from {current_location} to {self.pos}, and decreased the temperature there from {current_temperature:.2f} to {self.model.grid.properties['temperature'].data[self.pos]:.2f}"
        )

    def move_to_hottest_eligible_cell(self):
        # Define additional conditions
        conditions = {
            "humidity": lambda x: x > 30,  # Humidity in %
            "elevation": lambda x: x < 75,  # Elevation in meters
            "temperature": lambda x: x < 800,  # Temperature in Kelvin
        }

        # Create masks for the neighborhood
        neighborhood_mask = self.model.grid.get_neighborhood_mask(
            self.pos, moore=True, include_center=False, radius=3
        )

        # Select target cells based on conditions, extreme values, and masks
        target_cells = self.model.grid.select_cells(
            conditions=conditions,
            extreme_values={"temperature": "highest"},
            masks=neighborhood_mask,
            only_empty=True,
        )
        # Move the agent to one of the selected cells
        self.model.grid.move_agent_to_one_of(self, target_cells)

    def reduce_temperature(self):
        self.model.grid.properties["temperature"].modify_cell(
            self.pos, np.multiply, 0.8
        )


class ForestFireModel(Model):
    def __init__(self, N, width, height):
        super().__init__()
        self.num_agents = N
        self.grid = SingleGrid(width, height, False)
        self.schedule = RandomActivation(self)

        # Create property layers
        self.grid.add_property_layer(
            PropertyLayer("temperature", width, height, 300.0, dtype=float)
        )
        self.grid.add_property_layer(
            PropertyLayer("humidity", width, height, 60, dtype=int)
        )
        self.grid.add_property_layer(
            PropertyLayer("elevation", width, height, 0, dtype=int)
        )

        # Initialize property layers with random values
        self._init_property_layers()

        # Add agents
        for i in range(self.num_agents):
            agent = FirefighterAgent(i, self)
            self.grid.move_to_empty(agent)
            self.schedule.add(agent)

    def _init_property_layers(self):
        # Access the NumPy arrays of each property layer
        temperature_data = self.grid.properties["temperature"].data
        humidity_data = self.grid.properties["humidity"].data
        elevation_data = self.grid.properties["elevation"].data

        # Assign random values to the entire arrays at once
        temperature_data[:] = np.random.randint(270, 370, temperature_data.shape)
        humidity_data[:] = np.random.randint(20, 80, humidity_data.shape)
        elevation_data[:] = np.random.randint(0, 100, elevation_data.shape)

    def step(self):
        # Increase the temperature of all cells by 5
        print(
            f"Starting step {self.schedule.steps}. Average temperature: {self.grid.properties['temperature'].aggregate_property(np.mean):3f}"
        )
        self.grid.properties["temperature"].modify_cells(np.add, 5)

        self.schedule.step()


# Create and run the model
forest_fire_model = ForestFireModel(25, 100, 100)
for _i in range(100):  # Run for 10 steps
    forest_fire_model.step()
```

### Key Features of This Example

1. **Multiple Property Layers**: Temperature, humidity, and elevation layers model different aspects of the environment.
2. **Agent Behavior**: Firefighter agents move towards high-temperature zones while avoiding low humidity and high elevation areas, simulating strategic movement in a forest fire scenario.
3. **Dynamic Environment**: The model demonstrates how agents can interact with and respond to a multi-layered, dynamically changing environment.

codecov[bot]: ## [Codecov](https://app.codecov.io/gh/projectmesa/mesa/pull/1898?src=pr&el=h1&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=projectmesa) Report
Attention: `37 lines` in your changes are missing coverage. Please review.
> Comparison is base [(`838e216`)](https://app.codecov.io/gh/projectmesa/mesa/commit/838e21653d188ffbe54640d01b6a56fb3596244a?el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=projectmesa) 77.35% compared to head [(`221b4e3`)](https://app.codecov.io/gh/projectmesa/mesa/pull/1898?src=pr&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=projectmesa) 75.00%.

> :exclamation: Current head 221b4e3 differs from pull request most recent head e4be7d8. Consider uploading reports for the commit e4be7d8 to get more accurate results

| [Files](https://app.codecov.io/gh/projectmesa/mesa/pull/1898?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=projectmesa) | Patch % | Lines |
|---|---|---|
| [mesa/space.py](https://app.codecov.io/gh/projectmesa/mesa/pull/1898?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=projectmesa#diff-bWVzYS9zcGFjZS5weQ==) | 32.72% | [37 Missing :warning: ](https://app.codecov.io/gh/projectmesa/mesa/pull/1898?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=projectmesa) |

<details><summary>Additional details and impacted files</summary>


```diff
@@            Coverage Diff             @@
##             main    #1898      +/-   ##
==========================================
- Coverage   77.35%   75.00%   -2.36%     
==========================================
  Files          15       15              
  Lines        1007     1060      +53     
  Branches      220      230      +10     
==========================================
+ Hits          779      795      +16     
- Misses        197      234      +37     
  Partials       31       31              
```



</details>

[:umbrella: View full report in Codecov by Sentry](https://app.codecov.io/gh/projectmesa/mesa/pull/1898?src=pr&el=continue&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=projectmesa).   
:loudspeaker: Have feedback on the report? [Share it here](https://about.codecov.io/codecov-pr-comment-feedback/?utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=projectmesa).

EwoutH: This PR is now feature complete! Please take a proper look at the code, it's ready for review.

Next up:
- [ ] Tests
- [ ] Updating PR message
- [ ] Adding more examples
EwoutH: All docs and tests are added. That should reviewing a bit easier.

 - I'm quite certain I went for the right approach with PropertyLayer. I like the code base and the interface, and think it's quite performant and inuative
 - I also like the basic construct of adding one or more PropertyLayers to a grid. I think using a private _PropertyGrid of which SingleGrid and MultiGrid inherit is the right approach.
 - I'm not fully content with the build in methods of _PropertyGrid, especially the last three:
   - select_cells_multi_properties
   - move_agent_to_random_cell
   - move_agent_to_extreme_value_cell
   They work well but everything about them is just a little complicated, and maybe too specific. Any feedback on those will be extra appreciated.

While this PR is fully implemented and quite specific, feel also free to talk about the big picture stuff. I'm proposing a major conceptual change to the Mesa space stack, and thus discussion on all levels is appreciated.

<hr>

@rht Ruff doesn't allow me to specify a condition as a lambda function:

```python
condition = lambda x: x == 0
```
```bash
tests/test_space.py:557:9: E731 Do not assign a `lambda` expression, use a `def`
Found 7 errors.
No fixes available (7 hidden fixes can be enabled with the `--unsafe-fixes` option).
```
Direct insertion is okay, so without using the condition variable in between.

Being able to use lambda function to update the property layer is fundamental to this PR, and I like the current explicit `condition = lambda` approach. How do you think we should approach this. 
EwoutH: > One thing that got me thinking. You often define methods that potentially only affect a neighborhood. For this you always have the optional arguments only_neighborhoods, moore, include_center and radius. This somewhat predetermines what a neighborhood can be. But maybe we can define "neighborhood" just as a list of coordinates. Then we can define neighborhoods also based on conditions and only pass in the neighborhoods to these functions. Maybe that would be more general?

Yes, I like this. Just pass a sub selection of cells, which can be a neighbourhood but also can be something else. That would be extremely flexible, you can define rooms, areas, basically anything you like.
wang-boyu: Great work!

I'm wondering what type of indices is used for cells. In grid space the coordinate system is in `(x, y)` with origin at the lower left corner, which is typical Cartesian coordinates (if I remember correctly, but I might be wrong about this). However in matrices (or numpy arrays) we often use the `(row, column)` format, with an origin at the upper left corner. In mesa-geo I ended up using both: https://github.com/projectmesa/mesa-geo/blob/95cdcebb7490b0e832308ce8bb6161915c5de875/mesa_geo/raster_layers.py#L177C13-L180

Another question is that it seems to be possible to have property layers with different dimensions than that of the grid space. If this is the case, how do we align them? This also links to the coordinate system mentioned above. In mesa-geo, raster layers can have different dimensions and resolutions, and everything was taken care of through a `crs` setting in `GeoSpace`. But I think this is gis specific and may not apply here.

As an example, if you have a 10x10 grid, and a 20x20 property layer, where would you put it on the grid?

If, alternatively, you would prefer to make sure that property layers always have the same dimension as the grid, perhaps we could simplify the API to something like this?

```python
# Initialize a SingleGrid
grid = SingleGrid(width=10, height=10, torus=False)

# Create a PropertyLayer for a specific property (e.g., pollution level)
# and add the property layer to the grid
grid.add_property_layer(name='pollution', default_value=0.0)

# or
grid.add_property_layer(name='another_property', values=some_numpy_array_of_the_same_shape)
```
EwoutH: Based on Corvince’s insight and 
- Instead of neighborhood arguments, take a mask. It can be an empty cell mask, neighborhood mask, or any mask you like.
  - Add utility method to create mask from list of cells
  - Add utility methods to create neighborhood and empty cell masks
- Maybe split move_agent_to_extreme_value_cell into a select target and move to target method. That way you can mix and match them.
EwoutH: > I'm wondering what type of indices is used for cells. In grid space the coordinate system is in (x, y) with origin at the lower left corner, which is typical Cartesian coordinates (if I remember correctly, but I might be wrong about this). However in matrices (or numpy arrays) we often use the (row, column) format, with an origin at the upper left corner.

Not even thought of this. 😅

> Another question is that it seems to be possible to have property layers with different dimensions than that of the grid space. If this is the case, how do we align them?

Was also thinking about this for ContinuousSpace. I would say out of scope for this PR, but in general, it would be nice to support this.

> If, alternatively, you would prefer to make sure that property layers always have the same dimension as the grid, perhaps we could simplify the API to something like this?

Personally I wouldn’t go this way, because now property layers are completely independent. It would be nice to keep them that way.
EwoutH: Quite productive morning so far!

- In d731e2f I removed the neighborhood stuff from the PropertyGrid methods. They now take a mask argument
- In 478506b I allowed returning a mask or a list, for added flexibility
- In 5857c4f I split move_agent_to_extreme_value_cell into two methods: A select and a move method
- In c2f1b2d I wrote a utility function to create an empty mask
- In ca72bee I renamed two methods. I hope they're clearer now.
- In the other commits, I added docstring, updated the tests and fixed the remaining issues.

<hr>

@wang-boyu I added a test in 68992a0 to check the coordinate systems between agents and layers. It looks like they match. Can you check if I tested correctly and haven't forgotten an edge case?

<hr>

@rth I would really like to know what you think about Ruff [not allowing assigned lambda functions](https://github.com/projectmesa/mesa/pull/1898#issuecomment-1839508312).

<hr>

Code, test and docs wise, this PR is done. I will now work on some examples to show functionality. I would invite everyone to do that with me and play around with it a bit. This branch can be installed with:

```
pip install -U -e git+https://github.com/EwoutH/mesa@property_layer#egg=mesa
```
EwoutH: I added 3 examples to the main PR message. I think is shows (almost) all functionality. Really curious what everybody thinks!
wang-boyu: > @wang-boyu I added a test in [68992a0](https://github.com/projectmesa/mesa/commit/68992a0f3d634b45234932a7d9db3c9b467664f6) to check the coordinate systems between agents and layers. It looks like they match. Can you check if I tested correctly and haven't forgotten an edge case?

Thanks for adding a test case for it. But it's not really what I'm trying to say.

I made a small visualization for your forest fire model:

<details>
  <summary>code</summary>

```python
import numpy as np
import solara
import matplotlib.pyplot as plt
from matplotlib.figure import Figure
from mesa.experimental import JupyterViz


def space_drawer(model, agent_portrayal):
    agents_portrayal = {
        "firefighter": {"x": [], "y": [], "c": "tab:red", "marker": "o", "s": 10},
    }
    for agent in model.schedule.agents:
        if isinstance(agent, FirefighterAgent):
            agents_portrayal["firefighter"]["x"].append(agent.pos[0])
            agents_portrayal["firefighter"]["y"].append(agent.pos[1])

    fig, ax = plt.subplots()
    im = ax.imshow(model.grid.properties["temperature"].data, cmap="spring")
    # im = ax.imshow(model.grid.properties["temperature"].data, cmap="spring", origin="lower")
    fig.colorbar(im, orientation="vertical")

    ax.scatter(**agents_portrayal["firefighter"])
    ax.set_axis_off()
    solara.FigureMatplotlib(fig)


model_params = {
    "N": 1,
    "width": 4,
    "height": 4,
}

page = JupyterViz(
    ForestFireModel,
    model_params,
    measures=[],
    name="Forest Fire",
    space_drawer=space_drawer,
    play_interval=1500,
)
page  # noqa
```

</details>

For simplicity I also changed your `_init_property_layers()` function, from

```python
temperature_data[:] = np.random.randint(270, 370, temperature_data.shape)
```

to

```python
temperature_data[0, 0] = 500
```

so that the cell at position (0, 0) has a temperature of 500, with other cells having the default temperature of 300. Now the temperature property becomes

```python
model = ForestFireModel(1, 4, 4)
model.grid.properties["temperature"].data

array([[500., 300., 300., 300.],
       [300., 300., 300., 300.],
       [300., 300., 300., 300.],
       [300., 300., 300., 300.]])
```

If we plot it, it looks like

<img width="431" alt="Screenshot 2023-12-06 at 10 53 01" src="https://github.com/projectmesa/mesa/assets/10785873/1d1c60bd-469c-450e-9ac8-42a5d02e0cc1">

which isn't what we want, because (0, 0) is supposed to be at the lower left corner. This is essentially what I was trying to explain earlier.

I just realized that the same problem has been dealt with in the sugarscape demo: https://github.com/projectmesa/mesa-examples/blob/7275b172aad8905584dcffc131d9bbd96992abb7/examples/sugarscape_g1mt/app.py#L36-L38. In the demo code above I have a line commented out, that specifies `origin="lower"` when plotting the temperature layer. In that case, the visualization will look correct.

However I should mention that this only makes the model _looks_ correct. The underlying numpy array still has its _upper_ left corner with the 500 temperature. If we go for this approach, we probably need to make sure that

1. we and the users will use position in (x, y) format for consistency
2. all visualizations of property layers should have `origin="lower"`
3. the underlying numpy array is not the same as what it looks like in visualizations

Sorry about all these. I hope this helps clarifying my concern.
EwoutH: I already hate this, but adding a translation method that flips the y coordinate of every position like this:

```Python
    def _translate_coords(self, position):
        """
        Translates Mesa coordinates to NumPy array coordinates.
        """
        x, y = position
        return self.height - 1 - y, x
```
And then using that everywhere where a position is inputed:
```Python
    def set_cell(self, position, value):
        """
        Update a single cell's value in-place.
        """
        self.data[self._translate_coords(position)] = value
```
Would probably fix it right.

**Except** if users start manipulating the class directly, without using the methods...

I already hate this.
Corvince: Your hate is very understandable. What we could potentially do is use a numpy array with dtype=object for the agents as well. Last time I checked it was marginally slower, but would allow easier consistency. 
EwoutH: > We treat the numpy array as an implementation detail that is hidden from the users. Indexing is through (x, y) coordinate everywhere. We only need to make sure that, the right way to interpret the array (e.g., when visualizing it) is to set origin="lower". I believe this is what we did in the sugarscape example.

I thought about it for a bit, and I think this is the right approach. It's basically just a convention for visualisation. In future efforts, we might want to create a visualisation function that's 

> Regarding the first option, it is also possible to make the numpy array as a Python property without a setter.

I like this. Let me work on a implementation. We basically define it like this, right?

```Python
    @property
    def data(self):
        return self._data
```
<hr>

As a side note, I really would like feedback on the API and the variable names. This is the stuff that's difficult to change after the first release. So are there any names of methods that veel unintuitive, or could be improved?
EwoutH: To be honest, I don't know if we want data as a strict property. Doing some initializations as `np.arange(100).reshape(10, 10)`, `np.full((10, 10), 2)` is pretty powerful, and (as far as I understand) making `data` a property without a setter would disable that.

Or should we include a setter and do some checks there?
wang-boyu: > I thought about it for a bit, and I think this is the right approach. It's basically just a convention for visualisation. In future efforts, we might want to create a visualisation function that included the up-down flip.

Yes I prefer this approach too, which is also consistent with the sugarscape example model . To recap, we assume indexing such as `data[1, 2]` means `x=1` and `y=2`. This is different from the usual way of `(row, col)` indexing, but it is what we specify in Mesa.

A side effect of this approach is the internal representation using the numpy array will be upside down. To mitigate this problem, there are again a few options which are essentially design decisions:

1. Do nothing with the setter and getter. When plotting the layer, use `origin="lower"`. This is what we did in the sugarscape model.
2. Either transform the coordinate in the setter, or flip the numpy array in the getter. In this way the array representation will be correct, and we don't need to do `origin="lower"` in plots. An example of changing the getter would be (while the setter does nothing but returning the `_data`)
   ```python
   @property
   def data(self):
       return np.flipud(self._data)  # flip the layer; not very sure whether it works as expected
   ```
   So is the user does `data[0, 0] = 500`, then the cell at the _lower_ left corner gets updated. Might come as a surprise to those who are familiar with matrix indexing.

Regarding the APIs I don't have any preferences. My only questions is whether you would like to call it `cell` or `patch` or anything else.
Corvince: 
> Regarding the APIs I don't have any preferences. My only questions is whether you would like to call it `cell` or `patch` or anything else.

This is also something I noticed and I wonder if we should drop a concrete name altogether. I would be in favor of using more generic names for most of the functions. Anything with cell or even patch in the name implies a grid. If we later extend to 3 dimensional spaces or for networks the terminology doesn't fit anymore.

I would prefer to exchange cell with value in most cases. For example instead of `set_cell` we could have `set_value`, which fits for all kind of property layers. 

And `move_to_extreme_value_cell` could just drop the cell, which would also make it a bit shorter. 

In cases where value doesn't fit maybe we can use pos or position, which is still generic. 

Otherwise I am quite happy with this PR. Really excellent intro post, great examples.

Currently on mobile so I can't to a proper code review, but I agree that having a proper API is currently more important and implementation details can also be improved later (if there are any improvements)
EwoutH: Thanks! Let met think a bit about the `cell` / `patch` / `value` naming discussion. 

<hr>

As a demo, I've implemented Conway's Game of Life using the PropertyLayer:

```Python
import numpy as np
from mesa import Model
from mesa.space import SingleGrid, PropertyLayer
from scipy.signal import convolve2d

class GameOfLifeModel(Model):
    def __init__(self, width, height):
        self.grid = SingleGrid(width, height, torus=True)
        # Initialize the property layer for cell states
        self.cell_layer = PropertyLayer("cells", width, height, False, dtype=bool)
        self.grid.add_property_layer(self.cell_layer)

    from scipy.signal import convolve2d

    def step(self):
        # Define a kernel for counting neighbors. The kernel has 1s around the center cell (which is 0).
        # This setup allows us to count the live neighbors of each cell when we apply convolution.
        kernel = np.array([[1, 1, 1],
                           [1, 0, 1],
                           [1, 1, 1]])

        # Count neighbors using convolution.
        # convolve2d applies the kernel to each cell of the grid, summing up the values of neighbors.
        # boundary="wrap" ensures that the grid wraps around, simulating a toroidal surface.
        neighbor_count = convolve2d(self.cell_layer.data, kernel, mode="same", boundary="wrap")

        # Apply Game of Life rules:
        # 1. A live cell with 2 or 3 live neighbors survives, otherwise it dies.
        # 2. A dead cell with exactly 3 live neighbors becomes alive.
        # These rules are implemented using logical operations on the grid.
        self.cell_layer.data = np.logical_or(
            np.logical_and(self.cell_layer.data, np.logical_or(neighbor_count == 2, neighbor_count == 3)),
            # Rule for live cells
            np.logical_and(~self.cell_layer.data, neighbor_count == 3)  # Rule for dead cells
        )

    def count_live_neighbors(self, x, y):
        neighbors = self.grid.get_neighborhood((x, y), moore=True, include_center=False, radius=1)
        return sum(self.cell_layer.data[nx, ny] for nx, ny in neighbors)
```
EwoutH: I benchmarked it against [the GoL implementation](https://github.com/projectmesa/mesa-examples/tree/main/examples/conways_game_of_life) in the Mesa examples repo, and this implementation using a PropertyLayer is about 100x faster for the same grid size.

![runtime_comparison](https://github.com/projectmesa/mesa/assets/15776622/d30232c6-e23b-499b-8698-14695a95e627)

- [benchmark_gol.zip](https://github.com/projectmesa/mesa/files/13628343/benchmark_gol.zip)

EwoutH: I thought a bit about it, and while I like this discussion I have to disagree on this premise:

> Anything with cell or even patch in the name implies a grid.  If we later extend to 3 dimensional spaces

I think a cell implies a contained space, which can be 2D, 3D or otherwise. It also doesn’t imply a square grid, it could be hexagonal, multipolygon or else.

But it does imply a surface of volume is made up of discrete units, and each of the units takes some space, and all space is covered by exactly one unit.

> or for networks the terminology doesn't fit anymore.

I find networks difficult to be compatible with the PropertyLayer as implemented, because they don’t represent physical space. So I don’t think we’re ever going to mix those, since the PropertyLayer inherently implies physical space.

So personally, I currently prefer `cell`. I think it will also make teaching about it easier, there is something concrete you can talk about.
EwoutH: Initializing things randomly feels a bit weird:
```Python
class HabitatModel(Model):
    def __init__(self, n_prey, n_predators, food_probability):
        super().__init__()
        self.schedule = RandomActivation(self)
        self.running = True
        self.food_layer = PropertyLayer("food", 32, 32, default_value=False, dtype=bool)
        self.grid = MultiGrid(32, 32, True, property_layers=self.food_layer)
        self.food_probability = food_probability

        condition = lambda val: np.random.random_sample(val.shape) < self.food_probability
        self.food_layer.set_cells(True, condition)
```
I think random is important enough to built in. Doesn't have to be as part of this PR though.
wang-boyu: To recap my comments during today's dev meeting, which are mainly about how things are named. Apologies for the tediousness.

1. I'm also in favor of the following suggestion. Plus, `layer.data` could instead be `layer.values` to be consistent.

   > I would prefer to exchange cell with value in most cases. For example instead of `set_cell` we could have `set_value`, which fits for all kind of property layers.
   > 
   > And `move_to_extreme_value_cell` could just drop the cell, which would also make it a bit shorter.
   > 
   > In cases where value doesn't fit maybe we can use pos or position, which is still generic.

2. Regarding the function names:

   ```python
   set_cell(Coordinate position, value)
   set_cells(value, condition)
   ```
   Perhaps `set_cells()` could be `set_cells(condition, value)` (with parameters places swapped)? This is because `set_cell()` also has `value` as the second parameter.

   ```python
   modify_cell(Coordinate position, operation, value)
   modify_cells(operation, value, condition_function)
   ```

   This is where I miss function overloading in Python. I understand why these two functions exist. It's just that users need to remember when to use `set_cell()` and `modify_cell()` for similar things. Would it be too much to ask for something like `set_value(position, value=None, operation=None)`, and users are expected to provide either `value` or `operation`?

   For the case of both `operation` and `value` are needed in the current implementation, users can combine them as a single `operation`, maybe in the form of a lambda or partial function.

   We can even go further to just have one `set_value(position=None, condition=None, value=None, operation=None)` where one of `position` and `condition`, and one of `value` and `operation` is needed.

3. This PropertyLayer approach is much more efficient than before, and is equivalent to a simultaneous activation. We probably need to highlight this somewhere in the documentation, in case users aren't aware of it (e.g., if they really need a random activation for PropertyLayer, although I can't think of any particular example at the moment).
EwoutH: @quaquel curious what you think about this one. Going in the right direction, or shouldn't we entrench in (numpy) arrays to keep track of and manipulate environmental properties? Also considering future enhancements in data structures to the Grid itself.
quaquel: I haven't read through the entire discussion, so some points might have been discussed.

As I explained yesterday in person, I once rigged up my own SingleGrid implementation using NumPy. I had the grid itself and a 3d index array twice the size of the actual grid. Moore, von Neumann, and hex grids can all be done with this by changing how you slice on the indexing array.  Next, you flatten the return from the slice on the indexing array and have the neighborhood indices you then use on the actual grid to get the objects referenced here. I don't recall the details, but handling neighborhoods this way was several orders of magnitude faster. So, I am all in favor of building on NumPy for gridded spaces. 

I would avoid allowing property layers to have different gridding from the main grid. Likewise, I would standardize the origin to avoid the obnoxious translation between two different coordinate systems. By reimplementing grids themselves on top of NumPy, this problem would disappear.

Another thing to consider is that you might want to make it possible for PropertyLayers to have behavior. For example, to simulate grass regrowth. Conceptually, this would turn property layers into agents because they would have both state and behavior. So either property layers would subclass Agent, or you rely on duck typing through an Agent protocol to implement this cleanly. It also means you want to be able to schedule PropertyLayer updates. 


EwoutH: That pretty much aligns with the earlier discussion in this PR and in #1889. Only problem I see is even if you implement the Grid as a NumPy array, you still want to store a list of agents in each cell. Those data structured aren't optimized or vectorized, contrary to the numerical arrays here in the PropertyLayers.

> Another thing to consider is that you might want to make it possible for PropertyLayer to have behavior. 

In https://github.com/projectmesa/mesa/issues/1889#issuecomment-1831968467 we bisected this issues in two parts:
- A layer of patches that's passive and only has some property that can be read and modified.
- A proper Patch (or Cell) Agent

This PR solves the former, for the latter there isn't a proposal yet.

<hr>

As for this PR, I'm going to do a simple implementation of proposal #1903 to be able to move agents. Then we can remove the explicit movement functions from the `_PropertyGrid`.

@tpike3 and @jackiekazil after that I think this feature is in good enough state to be included in 2.2.0 as experimental.
quaquel: > That pretty much aligns with the earlier discussion in this PR and in #1889. Only problem I see is even if you implement the Grid as a NumPy array, you still want to store a list of agents in each cell. Those data structured aren't optimized or vectorized, contrary to the numerical arrays here in the PropertyLayers.
> 

True. However, the main performance difference comes from the slicing of the indexing array to get the indices of the neighboring cells. So you could still store stuff in the list of list things and only use NumPy to generate the indices given a neighborhood function, size, and whether you are on a torus or not. 

EwoutH: Updated move_agent method, which now can take a list of positions and choose either one randomly or the closest:

- #1920
EwoutH: With #1920 merged, I could remove the the move methods, making this PR a little less insanely big. 

```Python
target_cells = grid.select_cells_by_properties(conditions, mask=neighborhood_mask)
grid.move_agent_to_one_of(agent, target_cells)
```
Of course, you could make this an oneliner, but it might be a bit ugly.

```Python
grid.move_agent_to_one_of(agent, grid.select_cells_by_properties(conditions, mask=neighborhood_mask))
```

See the tests changed in a119cd9 for some more examples.

This PR is now fully ready for review, but please take your time, it's Christmas!🎄
EwoutH: @quaquel and @Corvince, any thoughts on my last two points?
quaquel: 1. sounds like a good idea
2. Not sure. Why take NetLogo as your source of inspiration? In essence a grid is a discretization, so why not try to support continuous properties were possible?
EwoutH: Will try an implementation for a single select function tomorrow.

> Not sure. Why take NetLogo as your source of inspiration?

We're standing on the shoulders of giants, right?

No mainly because many models and modelers are used to NetLogo, we should try to offer a superset of their capabilities. That makes porting models and enabling modelers to utilize Mesa easier.

> In essence a grid is a discretization, so why not try to support continuous properties were possible?

In the end that would be nice. Not sure how compatible that is with our current data structure (NumPy arrays). Maybe we can add some (optional) interpolation functionality when working with an continuous space.

I'm really trying to look hard for the the unknown unknowns.
- With a continuous space I imagine each agent having a (float, float) position. We will probably round (or floor) that when mapping to a cell as currently implemented in the PropertyLayer.
- When getting a value (so the other way around) we can either ask for the mapped value or an interpolated one.
- Setting values with some interpolation would probably be complicated, right?
- Radius checks (for neighbours etc.) could work like a float, and just check if the centers of patches are in it's radius

So I think the current PropertyLayer would work fine for most behavior. Still curious if I'm missing anything big.

We will be adding it as experimental, so maybe we should just start building it and then see if it works sufficiently or not.
Corvince: > With #1920 merged, I could remove the the move methods, making this PR a little less insanely big.
> 
> ```python
> target_cells = grid.select_cells_by_properties(conditions, mask=neighborhood_mask)
> grid.move_agent_to_one_of(agent, target_cells)
> ```
> 
> Of course, you could make this an oneliner, but it might be a bit ugly.
> 
> ```python
> grid.move_agent_to_one_of(agent, grid.select_cells_by_properties(conditions, mask=neighborhood_mask))
> ```
> 

Of course outside this PR, but this got me thinking. Maybe analog to our new AgentSet we can define something like a PositionSet and this could become 
```python
grid.select_cells_by_properties(conditions, mask=neighborhood_mask).move_to(agent)
```

Although even better to read would be 
´´´python
agent.move_to(grid.select_cells_by_properties(conditions, mask=neighborhood_mask).select(n=1))
```

But we haven't explored agents space function much yet. 

EwoutH: > Maybe analog to our new AgentSet we can define something like a PositionSet 

Also thought about it a few times. Each time really quickly tried to stop thinking about it due to the huge rewrite it would require.

But yeah, definitely worth thinking about.
EwoutH: I proposed in https://github.com/projectmesa/mesa/pull/1917#issuecomment-1870625103 to cut this feature out of the 2.2.0 release. We can discuss both the PositionSet as the ContiniousSpace implementation.
EwoutH: Thought about it some more, and I think it's best to first figure out how to port this to a `ContinousSpace`, and then think about a `PositionSet` kind of class. Because I think we first need to decide what a point is in a continuousspace and how to interpret that.

This means the HexGrid will be the only grid that doesn't work with the PropertyLayer. I don't like it, but I also don't see an logical way around it while keeping the PropertyLayer an ndarray.
quaquel: A set of points in continuous space is infinite by definition. So, if you want a `PositionSet` that works in continuous space, you have to somehow specify the boundary of the set. You cannot enumerate all positions, but in that way, you can at least check if points are in the set, and also generate any point that is within the set. I would suggest figuring out the API for the PositionSet first, starting for a simple grid. Once you have a fair idea of the API, we can worry about how to port this to a ContinuousSpace. 

What is the problem with hex spaces? In principle, you can implement hexspaces on top of an ndarray. you just get a slightly funny way of generating indices. I have done it a while ago when playing with Self Organizing Maps and learning numpy indexing and slicing. 


EwoutH: Time for a quick call?
EwoutH: Okay, I worked on this PR some more:
- 3d752a4 I merged the `select_cells_by_properties()` and `select_extreme_value_cells()` methods into a single, more flexible single `select_cells()` method that incorperates the functionality of both. It now also allows multiple extreme values and multiple masks. Further, it adds a boolean to only include empty cells. 
- 26d2c80 I replaced the slow `get_empty_mask()` method with a very fast `empty_mask` property. Also fixed a bug.
- 598bdea Marked the PropertyLayer as experimental. Feedback will be collected in #1932.

I also updated the UML and all the examples in the PR description.

I think this PR is ready for final review and to be merged. Let's see how it works in the field!

@jackiekazil @tpike3 @rht @Corvince @quaquel
EwoutH: Thanks for reviewing! I’m on vacation until the 7th, I will try to process your feedback the 8th.

Patch:
diff --git a/mesa/agent.py b/mesa/agent.py
index d62994ea..db452dbd 100644
--- a/mesa/agent.py
+++ b/mesa/agent.py
@@ -105,6 +105,8 @@ class AgentSet(MutableSet, Sequence):
         which means that agents not referenced elsewhere in the program may be automatically removed from the AgentSet.
     """
 
+    agentset_experimental_warning_given = False
+
     def __init__(self, agents: Iterable[Agent], model: Model):
         """
         Initializes the AgentSet with a collection of agents and a reference to the model.
@@ -115,8 +117,8 @@ class AgentSet(MutableSet, Sequence):
         """
         self.model = model
 
-        if not self.model.agentset_experimental_warning_given:
-            self.model.agentset_experimental_warning_given = True
+        if not self.__class__.agentset_experimental_warning_given:
+            self.__class__.agentset_experimental_warning_given = True
             warnings.warn(
                 "The AgentSet is experimental. It may be changed or removed in any and all future releases, including patch releases.\n"
                 "We would love to hear what you think about this new feature. If you have any thoughts, share them with us here: https://github.com/projectmesa/mesa/discussions/1919",
diff --git a/mesa/model.py b/mesa/model.py
index 6b5338da..ca0c653c 100644
--- a/mesa/model.py
+++ b/mesa/model.py
@@ -14,7 +14,7 @@ from collections import defaultdict
 # mypy
 from typing import Any
 
-from mesa.agent import AgentSet
+from mesa.agent import Agent, AgentSet
 from mesa.datacollection import DataCollector
 
 
@@ -83,7 +83,7 @@ class Model:
         """Return a list of different agent types."""
         return list(self._agents.keys())
 
-    def get_agents_of_type(self, agenttype: type) -> AgentSet:
+    def get_agents_of_type(self, agenttype: type[Agent]) -> AgentSet:
         """Retrieves an AgentSet containing all agents of the specified type."""
         return AgentSet(self._agents[agenttype].keys(), self)
 
diff --git a/mesa/space.py b/mesa/space.py
index 987afaf9..d1dc2ada 100644
--- a/mesa/space.py
+++ b/mesa/space.py
@@ -12,17 +12,16 @@ ContinuousSpace: a two-dimensional space where each agent has an arbitrary
                  position of `float`'s.
 NetworkGrid: a network where each node contains zero or more agents.
 """
-# Instruction for PyLint to suppress variable name errors, since we have a
-# good reason to use one-character variable names for x and y.
-# pylint: disable=invalid-name
 
 # Mypy; for the `|` operator purpose
 # Remove this __future__ import once the oldest supported Python is 3.10
 from __future__ import annotations
 
 import collections
+import inspect
 import itertools
 import math
+import warnings
 from numbers import Real
 from typing import (
     Any,
@@ -548,7 +547,424 @@ class _Grid:
         return len(self.empties) > 0
 
 
-class SingleGrid(_Grid):
+def is_single_argument_function(function):
+    """Check if a function is a single argument function."""
+    return (
+        inspect.isfunction(function)
+        and len(inspect.signature(function).parameters) == 1
+    )
+
+
+def ufunc_requires_additional_input(ufunc):
+    # NumPy ufuncs have a 'nargs' attribute indicating the number of input arguments
+    # For binary ufuncs (like np.add), nargs is 2
+    return ufunc.nargs > 1
+
+
+class PropertyLayer:
+    """
+    A class representing a layer of properties in a two-dimensional grid. Each cell in the grid
+    can store a value of a specified data type.
+
+    Attributes:
+        name (str): The name of the property layer.
+        width (int): The width of the grid (number of columns).
+        height (int): The height of the grid (number of rows).
+        data (numpy.ndarray): A NumPy array representing the grid data.
+
+    Methods:
+        set_cell(position, value): Sets the value of a single cell.
+        set_cells(value, condition=None): Sets the values of multiple cells, optionally based on a condition.
+        modify_cell(position, operation, value): Modifies the value of a single cell using an operation.
+        modify_cells(operation, value, condition_function): Modifies the values of multiple cells using an operation.
+        select_cells(condition, return_list): Selects cells that meet a specified condition.
+        aggregate_property(operation): Performs an aggregate operation over all cells.
+    """
+
+    agentset_experimental_warning_given = False
+
+    def __init__(
+        self, name: str, width: int, height: int, default_value, dtype=np.float64
+    ):
+        """
+        Initializes a new PropertyLayer instance.
+
+        Args:
+            name (str): The name of the property layer.
+            width (int): The width of the grid (number of columns). Must be a positive integer.
+            height (int): The height of the grid (number of rows). Must be a positive integer.
+            default_value: The default value to initialize each cell in the grid. Should ideally
+                           be of the same type as specified by the dtype parameter.
+            dtype (data-type, optional): The desired data-type for the grid's elements. Default is np.float64.
+
+        Raises:
+            ValueError: If width or height is not a positive integer.
+
+        Notes:
+            A UserWarning is raised if the default_value is not of a type compatible with dtype.
+            The dtype parameter can accept both Python data types (like bool, int or float) and NumPy data types
+            (like np.int64 or np.float64). Using NumPy data types is recommended (except for bool) for better control
+            over the precision and efficiency of data storage and computations, especially in cases of large data
+            volumes or specialized numerical operations.
+        """
+        self.name = name
+        self.width = width
+        self.height = height
+
+        # Check that width and height are positive integers
+        if (not isinstance(width, int) or width < 1) or (
+            not isinstance(height, int) or height < 1
+        ):
+            raise ValueError(
+                f"Width and height must be positive integers, got {width} and {height}."
+            )
+        # Check if the dtype is suitable for the data
+        if not isinstance(default_value, dtype):
+            warn(
+                f"Default value {default_value} ({type(default_value).__name__}) might not be best suitable with dtype={dtype.__name__}.",
+                UserWarning,
+                stacklevel=2,
+            )
+
+        self.data = np.full((width, height), default_value, dtype=dtype)
+
+        if not self.__class__.agentset_experimental_warning_given:
+            warnings.warn(
+                "The new PropertyLayer and _PropertyGrid classes experimental. It may be changed or removed in any and all future releases, including patch releases.\n"
+                "We would love to hear what you think about this new feature. If you have any thoughts, share them with us here: https://github.com/projectmesa/mesa/discussions/1932",
+                FutureWarning,
+                stacklevel=2,
+            )
+            self.__class__.agentset_experimental_warning_given = True
+
+    def set_cell(self, position: Coordinate, value):
+        """
+        Update a single cell's value in-place.
+        """
+        self.data[position] = value
+
+    def set_cells(self, value, condition=None):
+        """
+        Perform a batch update either on the entire grid or conditionally, in-place.
+
+        Args:
+            value: The value to be used for the update.
+            condition: (Optional) A callable (like a lambda function or a NumPy ufunc)
+                       that returns a boolean array when applied to the data.
+        """
+        if condition is None:
+            np.copyto(self.data, value)  # In-place update
+        else:
+            if isinstance(condition, np.ufunc):
+                # Directly apply NumPy ufunc
+                condition_result = condition(self.data)
+            else:
+                # Vectorize non-ufunc conditions
+                vectorized_condition = np.vectorize(condition)
+                condition_result = vectorized_condition(self.data)
+
+            if (
+                not isinstance(condition_result, np.ndarray)
+                or condition_result.shape != self.data.shape
+            ):
+                raise ValueError(
+                    "Result of condition must be a NumPy array with the same shape as the grid."
+                )
+
+            np.copyto(self.data, value, where=condition_result)
+
+    def modify_cell(self, position: Coordinate, operation, value=None):
+        """
+        Modify a single cell using an operation, which can be a lambda function or a NumPy ufunc.
+        If a NumPy ufunc is used, an additional value should be provided.
+
+        Args:
+            position: The grid coordinates of the cell to modify.
+            operation: A function to apply. Can be a lambda function or a NumPy ufunc.
+            value: The value to be used if the operation is a NumPy ufunc. Ignored for lambda functions.
+        """
+        current_value = self.data[position]
+
+        # Determine if the operation is a lambda function or a NumPy ufunc
+        if is_single_argument_function(operation):
+            # Lambda function case
+            self.data[position] = operation(current_value)
+        elif value is not None:
+            # NumPy ufunc case
+            self.data[position] = operation(current_value, value)
+        else:
+            raise ValueError("Invalid operation or missing value for NumPy ufunc.")
+
+    def modify_cells(self, operation, value=None, condition_function=None):
+        """
+        Modify cells using an operation, which can be a lambda function or a NumPy ufunc.
+        If a NumPy ufunc is used, an additional value should be provided.
+
+        Args:
+            operation: A function to apply. Can be a lambda function or a NumPy ufunc.
+            value: The value to be used if the operation is a NumPy ufunc. Ignored for lambda functions.
+            condition_function: (Optional) A callable that returns a boolean array when applied to the data.
+        """
+        condition_array = np.ones_like(
+            self.data, dtype=bool
+        )  # Default condition (all cells)
+        if condition_function is not None:
+            if isinstance(condition_function, np.ufunc):
+                condition_array = condition_function(self.data)
+            else:
+                vectorized_condition = np.vectorize(condition_function)
+                condition_array = vectorized_condition(self.data)
+
+        # Check if the operation is a lambda function or a NumPy ufunc
+        if isinstance(operation, np.ufunc):
+            if ufunc_requires_additional_input(operation):
+                if value is None:
+                    raise ValueError("This ufunc requires an additional input value.")
+                modified_data = operation(self.data, value)
+            else:
+                modified_data = operation(self.data)
+        else:
+            # Vectorize non-ufunc operations
+            vectorized_operation = np.vectorize(operation)
+            modified_data = vectorized_operation(self.data)
+
+        self.data = np.where(condition_array, modified_data, self.data)
+
+    def select_cells(self, condition, return_list=True):
+        """
+        Find cells that meet a specified condition using NumPy's boolean indexing, in-place.
+
+        Args:
+            condition: A callable that returns a boolean array when applied to the data.
+            return_list: (Optional) If True, return a list of (x, y) tuples. Otherwise, return a boolean array.
+
+        Returns:
+            A list of (x, y) tuples or a boolean array.
+        """
+        condition_array = condition(self.data)
+        if return_list:
+            return list(zip(*np.where(condition_array)))
+        else:
+            return condition_array
+
+    def aggregate_property(self, operation):
+        """Perform an aggregate operation (e.g., sum, mean) on a property across all cells.
+
+        Args:
+            operation: A function to apply. Can be a lambda function or a NumPy ufunc.
+        """
+        return operation(self.data)
+
+
+class _PropertyGrid(_Grid):
+    """
+    A private subclass of _Grid that supports the addition of property layers, enabling
+    the representation and manipulation of additional data layers on the grid. This class is
+    intended for internal use within the Mesa framework and is currently utilized by SingleGrid
+    and MultiGrid classes to provide enhanced grid functionality.
+
+    The `_PropertyGrid` extends the capabilities of a basic grid by allowing each cell
+    to have multiple properties, each represented by a separate PropertyLayer.
+    These properties can be used to model complex environments where each cell
+    has multiple attributes or states.
+
+    Attributes:
+        properties (dict): A dictionary mapping property layer names to PropertyLayer instances.
+        empty_mask: Returns a boolean mask indicating empty cells on the grid.
+
+    Methods:
+        add_property_layer(property_layer): Adds a new property layer to the grid.
+        remove_property_layer(property_name): Removes a property layer from the grid by its name.
+        get_neighborhood_mask(pos, moore, include_center, radius): Generates a boolean mask of the neighborhood.
+        select_cells_by_properties(conditions, mask, return_list): Selects cells based on multiple property conditions,
+            optionally with a mask, returning either a list of coordinates or a mask.
+        select_extreme_value_cells(property_name, mode, mask, return_list): Selects cells with extreme values of a property,
+            optionally with a mask, returning either a list of coordinates or a mask.
+
+    Mask Usage:
+        Several methods in this class accept a mask as an input, which is a NumPy ndarray of boolean values. This mask
+        specifies the cells to be considered (True) or ignored (False) in operations. Users can create custom masks,
+        including neighborhood masks, to apply specific conditions or constraints. Additionally, methods that deal with
+        cell selection or agent movement can return either a list of cell coordinates or a mask, based on the 'return_list'
+        parameter. This flexibility allows for more nuanced control and customization of grid operations, catering to a wide
+        range of modeling requirements and scenarios.
+
+    Note:
+        This class is not intended for direct use in user models but is currently used by the SingleGrid and MultiGrid.
+    """
+
+    def __init__(
+        self,
+        width: int,
+        height: int,
+        torus: bool,
+        property_layers: None | PropertyLayer | list[PropertyLayer] = None,
+    ):
+        """
+        Initializes a new _PropertyGrid instance with specified dimensions and optional property layers.
+
+        Args:
+            width (int): The width of the grid (number of columns).
+            height (int): The height of the grid (number of rows).
+            torus (bool): A boolean indicating if the grid should behave like a torus.
+            property_layers (None | PropertyLayer | list[PropertyLayer], optional): A single PropertyLayer instance,
+                a list of PropertyLayer instances, or None to initialize without any property layers.
+
+        Raises:
+            ValueError: If a property layer's dimensions do not match the grid dimensions.
+        """
+        super().__init__(width, height, torus)
+        self.properties = {}
+
+        # Initialize an empty mask as a boolean NumPy array
+        self._empty_mask = np.ones((self.width, self.height), dtype=bool)
+
+        # Handle both single PropertyLayer instance and list of PropertyLayer instances
+        if property_layers:
+            # If a single PropertyLayer is passed, convert it to a list
+            if isinstance(property_layers, PropertyLayer):
+                property_layers = [property_layers]
+
+            for layer in property_layers:
+                self.add_property_layer(layer)
+
+    @property
+    def empty_mask(self) -> np.ndarray:
+        """
+        Returns a boolean mask indicating empty cells on the grid.
+        """
+        return self._empty_mask
+
+    # Add and remove properties to the grid
+    def add_property_layer(self, property_layer: PropertyLayer):
+        """
+        Adds a new property layer to the grid.
+
+        Args:
+            property_layer (PropertyLayer): The PropertyLayer instance to be added to the grid.
+
+        Raises:
+            ValueError: If a property layer with the same name already exists in the grid.
+            ValueError: If the dimensions of the property layer do not match the grid's dimensions.
+        """
+        if property_layer.name in self.properties:
+            raise ValueError(f"Property layer {property_layer.name} already exists.")
+        if property_layer.width != self.width or property_layer.height != self.height:
+            raise ValueError(
+                f"Property layer dimensions {property_layer.width}x{property_layer.height} do not match grid dimensions {self.width}x{self.height}."
+            )
+        self.properties[property_layer.name] = property_layer
+
+    def remove_property_layer(self, property_name: str):
+        """
+        Removes a property layer from the grid by its name.
+
+        Args:
+            property_name (str): The name of the property layer to be removed.
+
+        Raises:
+            ValueError: If a property layer with the given name does not exist in the grid.
+        """
+        if property_name not in self.properties:
+            raise ValueError(f"Property layer {property_name} does not exist.")
+        del self.properties[property_name]
+
+    def get_neighborhood_mask(
+        self, pos: Coordinate, moore: bool, include_center: bool, radius: int
+    ) -> np.ndarray:
+        """
+        Generate a boolean mask representing the neighborhood.
+        Helper method for select_cells_multi_properties() and move_agent_to_random_cell()
+
+        Args:
+            pos (Coordinate): Center of the neighborhood.
+            moore (bool): True for Moore neighborhood, False for Von Neumann.
+            include_center (bool): Include the central cell in the neighborhood.
+            radius (int): The radius of the neighborhood.
+
+        Returns:
+            np.ndarray: A boolean mask representing the neighborhood.
+        """
+        neighborhood = self.get_neighborhood(pos, moore, include_center, radius)
+        mask = np.zeros((self.width, self.height), dtype=bool)
+
+        # Convert the neighborhood list to a NumPy array and use advanced indexing
+        coords = np.array(neighborhood)
+        mask[coords[:, 0], coords[:, 1]] = True
+        return mask
+
+    def select_cells(
+        self,
+        conditions: dict | None = None,
+        extreme_values: dict | None = None,
+        masks: np.ndarray | list[np.ndarray] = None,
+        only_empty: bool = False,
+        return_list: bool = True,
+    ) -> list[Coordinate] | np.ndarray:
+        """
+        Select cells based on property conditions, extreme values, and/or masks, with an option to only select empty cells.
+
+        Args:
+            conditions (dict): A dictionary where keys are property names and values are callables that return a boolean when applied.
+            extreme_values (dict): A dictionary where keys are property names and values are either 'highest' or 'lowest'.
+            masks (np.ndarray | list[np.ndarray], optional): A mask or list of masks to restrict the selection.
+            only_empty (bool, optional): If True, only select cells that are empty. Default is False.
+            return_list (bool, optional): If True, return a list of coordinates, otherwise return a mask.
+
+        Returns:
+            Union[list[Coordinate], np.ndarray]: Coordinates where conditions are satisfied or the combined mask.
+        """
+        # Initialize the combined mask
+        combined_mask = np.ones((self.width, self.height), dtype=bool)
+
+        # Apply the masks
+        if masks is not None:
+            if isinstance(masks, list):
+                for mask in masks:
+                    combined_mask = np.logical_and(combined_mask, mask)
+            else:
+                combined_mask = np.logical_and(combined_mask, masks)
+
+        # Apply the empty mask if only_empty is True
+        if only_empty:
+            combined_mask = np.logical_and(combined_mask, self.empty_mask)
+
+        # Apply conditions
+        if conditions:
+            for prop_name, condition in conditions.items():
+                prop_layer = self.properties[prop_name].data
+                prop_mask = condition(prop_layer)
+                combined_mask = np.logical_and(combined_mask, prop_mask)
+
+        # Apply extreme values
+        if extreme_values:
+            for property_name, mode in extreme_values.items():
+                prop_values = self.properties[property_name].data
+
+                # Create a masked array using the combined_mask
+                masked_values = np.ma.masked_array(prop_values, mask=~combined_mask)
+
+                if mode == "highest":
+                    target_value = masked_values.max()
+                elif mode == "lowest":
+                    target_value = masked_values.min()
+                else:
+                    raise ValueError(
+                        f"Invalid mode {mode}. Choose from 'highest' or 'lowest'."
+                    )
+
+                extreme_value_mask = prop_values == target_value
+                combined_mask = np.logical_and(combined_mask, extreme_value_mask)
+
+        # Generate output
+        if return_list:
+            selected_cells = list(zip(*np.where(combined_mask)))
+            return selected_cells
+        else:
+            return combined_mask
+
+
+class SingleGrid(_PropertyGrid):
     """Rectangular grid where each cell contains exactly at most one agent.
 
     Grid cells are indexed by [x, y], where [0, 0] is assumed to be the
@@ -576,6 +992,7 @@ class SingleGrid(_Grid):
             self._grid[x][y] = agent
             if self._empties_built:
                 self._empties.discard(pos)
+            self._empty_mask[pos] = False
             agent.pos = pos
         else:
             raise Exception("Cell not empty")
@@ -588,10 +1005,11 @@ class SingleGrid(_Grid):
         self._grid[x][y] = self.default_val()
         if self._empties_built:
             self._empties.add(pos)
+        self._empty_mask[agent.pos] = True
         agent.pos = None
 
 
-class MultiGrid(_Grid):
+class MultiGrid(_PropertyGrid):
     """Rectangular grid where each cell can contain more than one agent.
 
     Grid cells are indexed by [x, y], where [0, 0] is assumed to be at
@@ -623,6 +1041,7 @@ class MultiGrid(_Grid):
             agent.pos = pos
             if self._empties_built:
                 self._empties.discard(pos)
+                self._empty_mask[agent.pos] = True
 
     def remove_agent(self, agent: Agent) -> None:
         """Remove the agent from the given location and set its pos attribute to None."""
@@ -631,6 +1050,7 @@ class MultiGrid(_Grid):
         self._grid[x][y].remove(agent)
         if self._empties_built and self.is_cell_empty(pos):
             self._empties.add(pos)
+            self._empty_mask[agent.pos] = False
         agent.pos = None
 
     @accept_tuple_argument
diff --git a/mesa/time.py b/mesa/time.py
index c4f055eb..7bd4ad1b 100644
--- a/mesa/time.py
+++ b/mesa/time.py
@@ -27,12 +27,12 @@ from __future__ import annotations
 
 import heapq
 import warnings
-from collections import defaultdict
+import weakref
 
 # mypy
-from typing import Union
+from typing import Iterable, Union
 
-from mesa.agent import Agent
+from mesa.agent import Agent, AgentSet
 from mesa.model import Model
 
 # BaseScheduler has a self.time of int, while
@@ -60,12 +60,25 @@ class BaseScheduler:
         - agents (property): Returns a list of all agent instances.
     """
 
-    def __init__(self, model: Model) -> None:
-        """Create a new, empty BaseScheduler."""
+    def __init__(self, model: Model, agents: Iterable[Agent] | None = None) -> None:
+        """Create a new BaseScheduler.
+
+        Args:
+            model (Model): The model to which the schedule belongs
+            agents (Iterable[Agent], None, optional): An iterable of agents who are controlled by the schedule
+
+        """
         self.model = model
         self.steps = 0
         self.time: TimeT = 0
-        self._agents: dict[int, Agent] = {}
+
+        if agents is None:
+            agents = []
+
+        self._agents: AgentSet = AgentSet(agents, model)
+
+        self._remove_warning_given = False
+        self._agents_key_warning_given = False
 
     def add(self, agent: Agent) -> None:
         """Add an Agent object to the schedule.
@@ -74,20 +87,23 @@ class BaseScheduler:
             agent: An Agent to be added to the schedule. NOTE: The agent must
             have a step() method.
         """
-        if agent.unique_id in self._agents:
-            raise Exception(
-                f"Agent with unique id {agent.unique_id!r} already added to scheduler"
-            )
 
-        self._agents[agent.unique_id] = agent
+        if agent not in self._agents:
+            self._agents.add(agent)
+        else:
+            raise ValueError("agent already added to scheduler")
 
     def remove(self, agent: Agent) -> None:
         """Remove all instances of a given agent from the schedule.
 
+        Note:
+            It is only necessary to explicitly remove agents from the schedule if
+            the agent is not removed from the model.
+
         Args:
             agent: An agent object.
         """
-        del self._agents[agent.unique_id]
+        self._agents.remove(agent)
 
     def step(self) -> None:
         """Execute the step of all the agents, one at a time."""
@@ -102,25 +118,31 @@ class BaseScheduler:
         return len(self._agents)
 
     @property
-    def agents(self) -> list[Agent]:
-        return list(self._agents.values())
+    def agents(self) -> AgentSet:
+        # a bit dirty, but returns a copy of the internal agent set
+        return self._agents.select()
 
     def get_agent_keys(self, shuffle: bool = False) -> list[int]:
         # To be able to remove and/or add agents during stepping
         # it's necessary to cast the keys view to a list.
-        agent_keys = list(self._agents.keys())
+
+        if not self._agents_key_warning_given:
+            self._agents_key_warning_given = True
+            warnings.warn(
+                "Because of the shift to using weakrefs, this method will be removed in a future version",
+                DeprecationWarning,
+                stacklevel=2,
+            )
+
+        agent_keys = [agent.unique_id for agent in self._agents]
         if shuffle:
             self.model.random.shuffle(agent_keys)
         return agent_keys
 
-    def do_each(self, method, agent_keys=None, shuffle=False):
-        if agent_keys is None:
-            agent_keys = self.get_agent_keys()
+    def do_each(self, method, shuffle=False):
         if shuffle:
-            self.model.random.shuffle(agent_keys)
-        for agent_key in agent_keys:
-            if agent_key in self._agents:
-                getattr(self._agents[agent_key], method)()
+            self.agents.shuffle(inplace=True)
+        self.agents.do(method)
 
 
 class RandomActivation(BaseScheduler):
@@ -202,6 +224,7 @@ class StagedActivation(BaseScheduler):
     def __init__(
         self,
         model: Model,
+        agents: Iterable[Agent] | None = None,
         stage_list: list[str] | None = None,
         shuffle: bool = False,
         shuffle_between_stages: bool = False,
@@ -209,15 +232,16 @@ class StagedActivation(BaseScheduler):
         """Create an empty Staged Activation schedule.
 
         Args:
-            - model: Model object associated with the schedule.
-            - stage_list: List of strings of names of stages to run, in the
+            model (Model): The model to which the schedule belongs
+            agents (Iterable[Agent], None, optional): An iterable of agents who are controlled by the schedule
+            stage_list (:obj:`list` of :obj:`str`): List of strings of names of stages to run, in the
                          order to run them in.
-            - shuffle: If True, shuffle the order of agents each step.
-            - shuffle_between_stages: If True, shuffle the agents after each
+            shuffle (bool, optional): If True, shuffle the order of agents each step.
+            shuffle_between_stages (bool, optional): If True, shuffle the agents after each
                                     stage; otherwise, only shuffle at the start
                                     of each step.
         """
-        super().__init__(model)
+        super().__init__(model, agents)
         self.stage_list = stage_list if stage_list else ["step"]
         self.shuffle = shuffle
         self.shuffle_between_stages = shuffle_between_stages
@@ -225,17 +249,14 @@ class StagedActivation(BaseScheduler):
 
     def step(self) -> None:
         """Executes all the stages for all agents."""
-        # To be able to remove and/or add agents during stepping
-        # it's necessary for the keys view to be a list.
-        agent_keys = self.get_agent_keys(self.shuffle)
+        shuffle = self.shuffle
         for stage in self.stage_list:
             if stage.startswith("model."):
                 getattr(self.model, stage[6:])()
             else:
-                self.do_each(stage, agent_keys=agent_keys)
-            # We recompute the keys because some agents might have been removed
-            # in the previous loop.
-            agent_keys = self.get_agent_keys(self.shuffle_between_stages)
+                self.do_each(stage, shuffle=shuffle)
+
+            shuffle = self.shuffle_between_stages
             self.time += self.stage_time
 
         self.steps += 1
@@ -265,9 +286,24 @@ class RandomActivationByType(BaseScheduler):
         - get_type_count: Returns the count of agents of a specific type.
     """
 
-    def __init__(self, model: Model) -> None:
-        super().__init__(model)
-        self.agents_by_type = defaultdict(dict)
+    def __init__(self, model: Model, agents: Iterable[Agent] | None = None) -> None:
+        super().__init__(model, agents)
+        """
+
+        Args:
+            model (Model): The model to which the schedule belongs
+            agents (Iterable[Agent], None, optional): An iterable of agents who are controlled by the schedule
+        """
+
+        # can't be a defaultdict because we need to pass model to AgentSet
+        self.agents_by_type: [type, AgentSet] = {}
+
+        if agents is not None:
+            for agent in agents:
+                try:
+                    self.agents_by_type[type(agent)].add(agent)
+                except KeyError:
+                    self.agents_by_type[type(agent)] = AgentSet([agent], self.model)
 
     def add(self, agent: Agent) -> None:
         """
@@ -277,17 +313,18 @@ class RandomActivationByType(BaseScheduler):
             agent: An Agent to be added to the schedule.
         """
         super().add(agent)
-        agent_class: type[Agent] = type(agent)
-        self.agents_by_type[agent_class][agent.unique_id] = agent
+
+        try:
+            self.agents_by_type[type(agent)].add(agent)
+        except KeyError:
+            self.agents_by_type[type(agent)] = AgentSet([agent], self.model)
 
     def remove(self, agent: Agent) -> None:
         """
         Remove all instances of a given agent from the schedule.
         """
-        del self._agents[agent.unique_id]
-
-        agent_class: type[Agent] = type(agent)
-        del self.agents_by_type[agent_class][agent.unique_id]
+        super().remove(agent)
+        self.agents_by_type[type(agent)].remove(agent)
 
     def step(self, shuffle_types: bool = True, shuffle_agents: bool = True) -> None:
         """
@@ -309,26 +346,25 @@ class RandomActivationByType(BaseScheduler):
         self.steps += 1
         self.time += 1
 
-    def step_type(self, type_class: type[Agent], shuffle_agents: bool = True) -> None:
+    def step_type(self, agenttype: type[Agent], shuffle_agents: bool = True) -> None:
         """
         Shuffle order and run all agents of a given type.
         This method is equivalent to the NetLogo 'ask [breed]...'.
 
         Args:
-            type_class: Class object of the type to run.
+            agenttype: Class object of the type to run.
         """
-        agent_keys: list[int] = list(self.agents_by_type[type_class].keys())
+        agents = self.agents_by_type[agenttype]
+
         if shuffle_agents:
-            self.model.random.shuffle(agent_keys)
-        for agent_key in agent_keys:
-            if agent_key in self.agents_by_type[type_class]:
-                self.agents_by_type[type_class][agent_key].step()
+            agents.shuffle(inplace=True)
+        agents.do("step")
 
-    def get_type_count(self, type_class: type[Agent]) -> int:
+    def get_type_count(self, agenttype: type[Agent]) -> int:
         """
         Returns the current number of agents of certain type in the queue.
         """
-        return len(self.agents_by_type[type_class])
+        return len(self.agents_by_type[agenttype])
 
 
 class DiscreteEventScheduler(BaseScheduler):
@@ -376,8 +412,15 @@ class DiscreteEventScheduler(BaseScheduler):
     """
 
     def __init__(self, model: Model, time_step: TimeT = 1) -> None:
+        """
+
+        Args:
+            model (Model): The model to which the schedule belongs
+            time_step (TimeT): The fixed time step between steps
+
+        """
         super().__init__(model)
-        self.event_queue: list[tuple[TimeT, float, Agent]] = []
+        self.event_queue: list[tuple[TimeT, float, weakref.ref]] = []
         self.time_step: TimeT = time_step  # Fixed time period for each step
 
         warnings.warn(
@@ -393,8 +436,13 @@ class DiscreteEventScheduler(BaseScheduler):
             raise ValueError(
                 f"Scheduled time ({time}) must be >= the current time ({self.time})"
             )
+        if agent not in self._agents:
+            raise ValueError(
+                "trying to schedule an event for agent which is not known to the scheduler"
+            )
+
         # Create an event, sorted first on time, secondary on a random value
-        event = (time, self.model.random.random(), agent)
+        event = (time, self.model.random.random(), weakref.ref(agent))
         heapq.heappush(self.event_queue, event)
 
     def schedule_in(self, delay: TimeT, agent: Agent) -> None:
@@ -411,12 +459,12 @@ class DiscreteEventScheduler(BaseScheduler):
         while self.event_queue and self.event_queue[0][0] <= end_time:
             # Get the next event (ignore the random value during unpacking)
             time, _, agent = heapq.heappop(self.event_queue)
+            agent = agent()  # unpack weakref
 
-            # Advance time to the event's time
-            self.time = time
-
-            # Execute the event
-            if agent.unique_id in self._agents:
+            if agent:
+                # Advance time to the event's time
+                self.time = time
+                # Execute the event
                 agent.step()
 
         # After processing events, advance time by the time_step


Test patch:
diff --git a/tests/test_space.py b/tests/test_space.py
index ab0e0da3..b7524b6d 100644
--- a/tests/test_space.py
+++ b/tests/test_space.py
@@ -4,7 +4,7 @@ import networkx as nx
 import numpy as np
 import pytest
 
-from mesa.space import ContinuousSpace, NetworkGrid, SingleGrid
+from mesa.space import ContinuousSpace, NetworkGrid, PropertyLayer, SingleGrid
 from tests.test_grid import MockAgent
 
 TEST_AGENTS = [(-20, -20), (-20, -20.05), (65, 18)]
@@ -278,6 +278,165 @@ class TestSpaceAgentMapping(unittest.TestCase):
             self.space.remove_agent(agent_to_remove)
 
 
+class TestPropertyLayer(unittest.TestCase):
+    def setUp(self):
+        self.layer = PropertyLayer("test_layer", 10, 10, 0)
+
+    # Initialization Test
+    def test_initialization(self):
+        self.assertEqual(self.layer.name, "test_layer")
+        self.assertEqual(self.layer.width, 10)
+        self.assertEqual(self.layer.height, 10)
+        self.assertTrue(np.array_equal(self.layer.data, np.zeros((10, 10))))
+
+    # Set Cell Test
+    def test_set_cell(self):
+        self.layer.set_cell((5, 5), 1)
+        self.assertEqual(self.layer.data[5, 5], 1)
+
+    # Set Cells Tests
+    def test_set_cells_no_condition(self):
+        self.layer.set_cells(2)
+        np.testing.assert_array_equal(self.layer.data, np.full((10, 10), 2))
+
+    def test_set_cells_with_condition(self):
+        self.layer.set_cell((5, 5), 1)
+
+        def condition(x):
+            return x == 0
+
+        self.layer.set_cells(3, condition)
+        self.assertEqual(self.layer.data[5, 5], 1)
+        self.assertEqual(self.layer.data[0, 0], 3)
+        # Check if the sum is correct
+        self.assertEqual(np.sum(self.layer.data), 3 * 99 + 1)
+
+    def test_set_cells_with_random_condition(self):
+        # Probability for a cell to be updated
+        update_probability = 0.5
+
+        # Define a condition with a random part
+        def condition(val):
+            return np.random.rand() < update_probability
+
+        # Apply set_cells
+        self.layer.set_cells(True, condition)
+
+        # Count the number of cells that were set to True
+        true_count = np.sum(self.layer.data)
+
+        width = self.layer.width
+        height = self.layer.height
+
+        # Calculate expected range (with some tolerance for randomness)
+        expected_min = width * height * update_probability * 0.5
+        expected_max = width * height * update_probability * 1.5
+
+        # Check if the true_count falls within the expected range
+        assert expected_min <= true_count <= expected_max
+
+    # Modify Cell Test
+    def test_modify_cell_lambda(self):
+        self.layer.data = np.zeros((10, 10))
+        self.layer.modify_cell((2, 2), lambda x: x + 5)
+        self.assertEqual(self.layer.data[2, 2], 5)
+
+    def test_modify_cell_ufunc(self):
+        self.layer.data = np.ones((10, 10))
+        self.layer.modify_cell((3, 3), np.add, 4)
+        self.assertEqual(self.layer.data[3, 3], 5)
+
+    def test_modify_cell_invalid_operation(self):
+        with self.assertRaises(ValueError):
+            self.layer.modify_cell((1, 1), np.add)  # Missing value for ufunc
+
+    # Modify Cells Test
+    def test_modify_cells_lambda(self):
+        self.layer.data = np.zeros((10, 10))
+        self.layer.modify_cells(lambda x: x + 2)
+        np.testing.assert_array_equal(self.layer.data, np.full((10, 10), 2))
+
+    def test_modify_cells_ufunc(self):
+        self.layer.data = np.ones((10, 10))
+        self.layer.modify_cells(np.multiply, 3)
+        np.testing.assert_array_equal(self.layer.data, np.full((10, 10), 3))
+
+    def test_modify_cells_invalid_operation(self):
+        with self.assertRaises(ValueError):
+            self.layer.modify_cells(np.add)  # Missing value for ufunc
+
+    # Aggregate Property Test
+    def test_aggregate_property_lambda(self):
+        self.layer.data = np.arange(100).reshape(10, 10)
+        result = self.layer.aggregate_property(lambda x: np.sum(x))
+        self.assertEqual(result, np.sum(np.arange(100)))
+
+    def test_aggregate_property_ufunc(self):
+        self.layer.data = np.full((10, 10), 2)
+        result = self.layer.aggregate_property(np.mean)
+        self.assertEqual(result, 2)
+
+    # Edge Case: Negative or Zero Dimensions
+    def test_initialization_negative_dimensions(self):
+        with self.assertRaises(ValueError):
+            PropertyLayer("test_layer", -10, 10, 0)
+
+    def test_initialization_zero_dimensions(self):
+        with self.assertRaises(ValueError):
+            PropertyLayer("test_layer", 0, 10, 0)
+
+    # Edge Case: Out-of-Bounds Cell Access
+    def test_set_cell_out_of_bounds(self):
+        with self.assertRaises(IndexError):
+            self.layer.set_cell((10, 10), 1)
+
+    def test_modify_cell_out_of_bounds(self):
+        with self.assertRaises(IndexError):
+            self.layer.modify_cell((10, 10), lambda x: x + 5)
+
+    # Edge Case: Selecting Cells with Complex Conditions
+    def test_select_cells_complex_condition(self):
+        self.layer.data = np.random.rand(10, 10)
+        selected = self.layer.select_cells(lambda x: (x > 0.5) & (x < 0.75))
+        for c in selected:
+            self.assertTrue(0.5 < self.layer.data[c] < 0.75)
+
+    # More edge cases
+    def test_set_cells_with_numpy_ufunc(self):
+        # Set some cells to a specific value
+        self.layer.data[0:5, 0:5] = 5
+
+        # Use a numpy ufunc as a condition. Here, we're using `np.greater`
+        # which will return True for cells with values greater than 2.
+        condition = np.greater
+        self.layer.set_cells(10, lambda x: condition(x, 2))
+
+        # Check if cells that had value greater than 2 are now set to 10
+        updated_cells = self.layer.data[0:5, 0:5]
+        np.testing.assert_array_equal(updated_cells, np.full((5, 5), 10))
+
+        # Check if cells that had value 0 (less than or equal to 2) remain unchanged
+        unchanged_cells = self.layer.data[5:, 5:]
+        np.testing.assert_array_equal(unchanged_cells, np.zeros((5, 5)))
+
+    def test_modify_cell_boundary_condition(self):
+        self.layer.data = np.zeros((10, 10))
+        self.layer.modify_cell((0, 0), lambda x: x + 5)
+        self.layer.modify_cell((9, 9), lambda x: x + 5)
+        self.assertEqual(self.layer.data[0, 0], 5)
+        self.assertEqual(self.layer.data[9, 9], 5)
+
+    def test_aggregate_property_std_dev(self):
+        self.layer.data = np.arange(100).reshape(10, 10)
+        result = self.layer.aggregate_property(np.std)
+        self.assertAlmostEqual(result, np.std(np.arange(100)), places=5)
+
+    def test_data_type_consistency(self):
+        self.layer.data = np.zeros((10, 10), dtype=int)
+        self.layer.set_cell((5, 5), 5.5)
+        self.assertIsInstance(self.layer.data[5, 5], self.layer.data.dtype.type)
+
+
 class TestSingleGrid(unittest.TestCase):
     def setUp(self):
         self.space = SingleGrid(50, 50, False)
@@ -312,6 +471,16 @@ class TestSingleGrid(unittest.TestCase):
         with self.assertRaises(Exception):
             self.space.move_to_empty(a)
 
+    def test_empty_mask_consistency(self):
+        # Check that the empty mask is consistent with the empties set
+        empty_mask = self.space.empty_mask
+        empties = self.space.empties
+        for i in range(self.space.width):
+            for j in range(self.space.height):
+                mask_value = empty_mask[i, j]
+                empties_value = (i, j) in empties
+                assert mask_value == empties_value
+
     def move_agent(self):
         agent_number = 0
         initial_pos = TEST_AGENTS_GRID[agent_number]
@@ -439,6 +608,242 @@ class TestSingleGridTorus(unittest.TestCase):
         assert self.space._distance_squared(pos1, pos2) == expected_distance_squared
 
 
+class TestSingleGridWithPropertyGrid(unittest.TestCase):
+    def setUp(self):
+        self.grid = SingleGrid(10, 10, False)
+        self.property_layer1 = PropertyLayer("layer1", 10, 10, 0)
+        self.property_layer2 = PropertyLayer("layer2", 10, 10, 1)
+        self.grid.add_property_layer(self.property_layer1)
+        self.grid.add_property_layer(self.property_layer2)
+
+    # Test adding and removing property layers
+    def test_add_property_layer(self):
+        self.assertIn("layer1", self.grid.properties)
+        self.assertIn("layer2", self.grid.properties)
+
+    def test_remove_property_layer(self):
+        self.grid.remove_property_layer("layer1")
+        self.assertNotIn("layer1", self.grid.properties)
+
+    def test_add_property_layer_mismatched_dimensions(self):
+        with self.assertRaises(ValueError):
+            self.grid.add_property_layer(PropertyLayer("layer3", 5, 5, 0))
+
+    def test_add_existing_property_layer(self):
+        with self.assertRaises(ValueError):
+            self.grid.add_property_layer(self.property_layer1)
+
+    def test_remove_nonexistent_property_layer(self):
+        with self.assertRaises(ValueError):
+            self.grid.remove_property_layer("nonexistent_layer")
+
+    # Test getting masks
+    def test_get_empty_mask(self):
+        empty_mask = self.grid.empty_mask
+        self.assertTrue(np.all(empty_mask == np.ones((10, 10), dtype=bool)))
+
+    def test_get_empty_mask_with_agent(self):
+        agent = MockAgent(0, self.grid)
+        self.grid.place_agent(agent, (4, 6))
+
+        empty_mask = self.grid.empty_mask
+        expected_mask = np.ones((10, 10), dtype=bool)
+        expected_mask[4, 6] = False
+
+        self.assertTrue(np.all(empty_mask == expected_mask))
+
+    def test_get_neighborhood_mask(self):
+        agent = MockAgent(0, self.grid)
+        agent2 = MockAgent(1, self.grid)
+        self.grid.place_agent(agent, (5, 5))
+        self.grid.place_agent(agent2, (5, 6))
+        neighborhood_mask = self.grid.get_neighborhood_mask((5, 5), True, False, 1)
+        expected_mask = np.zeros((10, 10), dtype=bool)
+        expected_mask[4:7, 4:7] = True
+        expected_mask[5, 5] = False
+        self.assertTrue(np.all(neighborhood_mask == expected_mask))
+
+    # Test selecting and moving to cells based on multiple conditions
+    def test_select_cells_by_properties(self):
+        def condition(x):
+            return x == 0
+
+        selected_cells = self.grid.select_cells({"layer1": condition})
+        self.assertEqual(len(selected_cells), 100)
+
+    def test_select_cells_by_properties_return_mask(self):
+        def condition(x):
+            return x == 0
+
+        selected_mask = self.grid.select_cells({"layer1": condition}, return_list=False)
+        self.assertTrue(isinstance(selected_mask, np.ndarray))
+        self.assertTrue(selected_mask.all())
+
+    def test_move_agent_to_cell_by_properties(self):
+        agent = MockAgent(1, self.grid)
+        self.grid.place_agent(agent, (5, 5))
+        conditions = {"layer1": lambda x: x == 0}
+        target_cells = self.grid.select_cells(conditions)
+        self.grid.move_agent_to_one_of(agent, target_cells)
+        # Agent should move, since none of the cells match the condition
+        self.assertNotEqual(agent.pos, (5, 5))
+
+    def test_move_agent_no_eligible_cells(self):
+        agent = MockAgent(3, self.grid)
+        self.grid.place_agent(agent, (5, 5))
+        conditions = {"layer1": lambda x: x != 0}
+        target_cells = self.grid.select_cells(conditions)
+        self.grid.move_agent_to_one_of(agent, target_cells)
+        self.assertEqual(agent.pos, (5, 5))
+
+    # Test selecting and moving to cells based on extreme values
+    def test_select_extreme_value_cells(self):
+        self.grid.properties["layer2"].set_cell((3, 1), 1.1)
+        target_cells = self.grid.select_cells(extreme_values={"layer2": "highest"})
+        self.assertIn((3, 1), target_cells)
+
+    def test_select_extreme_value_cells_return_mask(self):
+        self.grid.properties["layer2"].set_cell((3, 1), 1.1)
+        target_mask = self.grid.select_cells(
+            extreme_values={"layer2": "highest"}, return_list=False
+        )
+        self.assertTrue(isinstance(target_mask, np.ndarray))
+        self.assertTrue(target_mask[3, 1])
+
+    def test_move_agent_to_extreme_value_cell(self):
+        agent = MockAgent(2, self.grid)
+        self.grid.place_agent(agent, (5, 5))
+        self.grid.properties["layer2"].set_cell((3, 1), 1.1)
+        target_cells = self.grid.select_cells(extreme_values={"layer2": "highest"})
+        self.grid.move_agent_to_one_of(agent, target_cells)
+        self.assertEqual(agent.pos, (3, 1))
+
+    # Test using masks
+    def test_select_cells_by_properties_with_empty_mask(self):
+        self.grid.place_agent(
+            MockAgent(0, self.grid), (5, 5)
+        )  # Placing an agent to ensure some cells are not empty
+        empty_mask = self.grid.empty_mask
+
+        def condition(x):
+            return x == 0
+
+        selected_cells = self.grid.select_cells({"layer1": condition}, masks=empty_mask)
+        self.assertNotIn(
+            (5, 5), selected_cells
+        )  # (5, 5) should not be in the selection as it's not empty
+
+    def test_select_cells_by_properties_with_neighborhood_mask(self):
+        neighborhood_mask = self.grid.get_neighborhood_mask((5, 5), True, False, 1)
+
+        def condition(x):
+            return x == 0
+
+        selected_cells = self.grid.select_cells(
+            {"layer1": condition}, masks=neighborhood_mask
+        )
+        expected_selection = [
+            (4, 4),
+            (4, 5),
+            (4, 6),
+            (5, 4),
+            (5, 6),
+            (6, 4),
+            (6, 5),
+            (6, 6),
+        ]  # Cells in the neighborhood of (5, 5)
+        self.assertCountEqual(selected_cells, expected_selection)
+
+    def test_move_agent_to_cell_by_properties_with_empty_mask(self):
+        agent = MockAgent(1, self.grid)
+        self.grid.place_agent(agent, (5, 5))
+        self.grid.place_agent(
+            MockAgent(2, self.grid), (4, 5)
+        )  # Placing another agent to create a non-empty cell
+        empty_mask = self.grid.empty_mask
+        conditions = {"layer1": lambda x: x == 0}
+        target_cells = self.grid.select_cells(conditions, masks=empty_mask)
+        self.grid.move_agent_to_one_of(agent, target_cells)
+        self.assertNotEqual(
+            agent.pos, (4, 5)
+        )  # Agent should not move to (4, 5) as it's not empty
+
+    def test_move_agent_to_cell_by_properties_with_neighborhood_mask(self):
+        agent = MockAgent(1, self.grid)
+        self.grid.place_agent(agent, (5, 5))
+        neighborhood_mask = self.grid.get_neighborhood_mask((5, 5), True, False, 1)
+        conditions = {"layer1": lambda x: x == 0}
+        target_cells = self.grid.select_cells(conditions, masks=neighborhood_mask)
+        self.grid.move_agent_to_one_of(agent, target_cells)
+        self.assertIn(
+            agent.pos, [(4, 4), (4, 5), (4, 6), (5, 4), (5, 6), (6, 4), (6, 5), (6, 6)]
+        )  # Agent should move within the neighborhood
+
+    # Test invalid inputs
+    def test_invalid_property_name_in_conditions(self):
+        def condition(x):
+            return x == 0
+
+        with self.assertRaises(KeyError):
+            self.grid.select_cells(conditions={"nonexistent_layer": condition})
+
+    # Test if coordinates means the same between the grid and the property layer
+    def test_property_layer_coordinates(self):
+        agent = MockAgent(0, self.grid)
+        correct_pos = (1, 8)
+        incorrect_pos = (8, 1)
+        self.grid.place_agent(agent, correct_pos)
+
+        # Simple check on layer 1: set by agent, check by layer
+        self.grid.properties["layer1"].set_cell(agent.pos, 2)
+        self.assertEqual(self.grid.properties["layer1"].data[agent.pos], 2)
+
+        # More complicated check on layer 2: set by layer, check by agent
+        self.grid.properties["layer2"].set_cell(correct_pos, 3)
+        self.grid.properties["layer2"].set_cell(incorrect_pos, 4)
+
+        correct_grid_value = self.grid.properties["layer2"].data[correct_pos]
+        incorrect_grid_value = self.grid.properties["layer2"].data[incorrect_pos]
+        agent_grid_value = self.grid.properties["layer2"].data[agent.pos]
+
+        self.assertEqual(correct_grid_value, agent_grid_value)
+        self.assertNotEqual(incorrect_grid_value, agent_grid_value)
+
+    # Test selecting cells with only_empty parameter
+    def test_select_cells_only_empty(self):
+        self.grid.place_agent(MockAgent(0, self.grid), (5, 5))  # Occupying a cell
+        selected_cells = self.grid.select_cells(only_empty=True)
+        self.assertNotIn(
+            (5, 5), selected_cells
+        )  # The occupied cell should not be selected
+
+    def test_select_cells_only_empty_with_conditions(self):
+        self.grid.place_agent(MockAgent(1, self.grid), (5, 5))
+        self.grid.properties["layer1"].set_cell((5, 5), 2)
+        self.grid.properties["layer1"].set_cell((6, 6), 2)
+
+        def condition(x):
+            return x == 2
+
+        selected_cells = self.grid.select_cells({"layer1": condition}, only_empty=True)
+        self.assertIn((6, 6), selected_cells)
+        self.assertNotIn((5, 5), selected_cells)
+
+    # Test selecting cells with multiple extreme values
+    def test_select_cells_multiple_extreme_values(self):
+        self.grid.properties["layer1"].set_cell((1, 1), 3)
+        self.grid.properties["layer1"].set_cell((2, 2), 3)
+        self.grid.properties["layer2"].set_cell((2, 2), 0.5)
+        self.grid.properties["layer2"].set_cell((3, 3), 0.5)
+        selected_cells = self.grid.select_cells(
+            extreme_values={"layer1": "highest", "layer2": "lowest"}
+        )
+        self.assertIn((2, 2), selected_cells)
+        self.assertNotIn((1, 1), selected_cells)
+        self.assertNotIn((3, 3), selected_cells)
+        self.assertEqual(len(selected_cells), 1)
+
+
 class TestSingleNetworkGrid(unittest.TestCase):
     GRAPH_SIZE = 10
 
diff --git a/tests/test_time.py b/tests/test_time.py
index 5d6d66f9..a440d502 100644
--- a/tests/test_time.py
+++ b/tests/test_time.py
@@ -34,8 +34,7 @@ class MockAgent(Agent):
     def kill_other_agent(self):
         for agent in self.model.schedule.agents:
             if agent is not self:
-                self.model.schedule.remove(agent)
-                break
+                agent.remove()
 
     def stage_one(self):
         if self.model.enable_kill_other_agent:
@@ -78,7 +77,9 @@ class MockModel(Model):
         # Make scheduler
         if activation == STAGED:
             model_stages = ["stage_one", "model.model_stage", "stage_two"]
-            self.schedule = StagedActivation(self, model_stages, shuffle=shuffle)
+            self.schedule = StagedActivation(
+                self, stage_list=model_stages, shuffle=shuffle
+            )
         elif activation == RANDOM:
             self.schedule = RandomActivation(self)
         elif activation == SIMULTANEOUS:
@@ -140,9 +141,9 @@ class TestStagedActivation(TestCase):
         Test the staged activation can remove an agent
         """
         model = MockModel(shuffle=True)
-        agent_keys = list(model.schedule._agents.keys())
-        agent = model.schedule._agents[agent_keys[0]]
-        model.schedule.remove(agent)
+        agents = list(model.schedule._agents)
+        agent = agents[0]
+        model.schedule.remove(agents[0])
         assert agent not in model.schedule.agents
 
     def test_intrastep_remove(self):
@@ -166,6 +167,14 @@ class TestRandomActivation(TestCase):
     Test the random activation.
     """
 
+    def test_init(self):
+        model = Model()
+        agents = [MockAgent(model.next_id(), model) for _ in range(10)]
+
+        scheduler = RandomActivation(model, agents)
+        assert all(agent in scheduler.agents for agent in agents)
+
+
     def test_random_activation_step_shuffles(self):
         """
         Test the random activation step
@@ -205,6 +214,18 @@ class TestRandomActivation(TestCase):
         model.step()
         assert len(model.log) == 1
 
+    def test_get_agent_keys(self):
+        model = MockModel(activation=RANDOM)
+
+        keys = model.schedule.get_agent_keys()
+        agent_ids = [agent.unique_id for agent in model.agents]
+        assert all(entry_i == entry_j for entry_i, entry_j in zip(keys, agent_ids))
+
+        keys = model.schedule.get_agent_keys(shuffle=True)
+        agent_ids = {agent.unique_id for agent in model.agents}
+        assert all(entry in agent_ids for entry in keys)
+
+
 
 class TestSimultaneousActivation(TestCase):
     """
@@ -231,6 +252,15 @@ class TestRandomActivationByType(TestCase):
     does step for one type of agents, not the entire agents.
     """
 
+    def test_init(self):
+        model = Model()
+        agents = [MockAgent(model.next_id(), model) for _ in range(10)]
+        agents += [Agent(model.next_id(), model) for _ in range(10)]
+
+        scheduler = RandomActivationByType(model, agents)
+        assert all(agent in scheduler.agents for agent in agents)
+
+
     def test_random_activation_step_shuffles(self):
         """
         Test the random activation by type step
@@ -262,18 +292,30 @@ class TestRandomActivationByType(TestCase):
         # one step for each of 2 agents
         assert all(x == 1 for x in agent_steps)
 
-    def test_add_non_unique_ids(self):
+    def test_random_activation_counts(self):
         """
-        Test that adding agent with duplicate ids result in an error.
-        TODO: we need to run this test on all schedulers, not just
-        RandomActivationByType.
+        Test the random activation by type step causes each agent to step
         """
+
         model = MockModel(activation=RANDOM_BY_TYPE)
-        a = MockAgent(0, model)
-        b = MockAgent(0, model)
-        model.schedule.add(a)
-        with self.assertRaises(Exception):
-            model.schedule.add(b)
+
+        agent_types = model.agent_types
+        for agent_type in agent_types:
+            assert model.schedule.get_type_count(agent_type) == len(model.get_agents_of_type(agent_type))
+
+    # def test_add_non_unique_ids(self):
+    #     """
+    #     Test that adding agent with duplicate ids result in an error.
+    #     TODO: we need to run this test on all schedulers, not just
+    #     TODO:: identical IDs is something for the agent, not the scheduler and should be tested there
+    #     RandomActivationByType.
+    #     """
+    #     model = MockModel(activation=RANDOM_BY_TYPE)
+    #     a = MockAgent(0, model)
+    #     b = MockAgent(0, model)
+    #     model.schedule.add(a)
+    #     with self.assertRaises(Exception):
+    #         model.schedule.add(b)
 
 
 class TestDiscreteEventScheduler(TestCase):
@@ -297,19 +339,19 @@ class TestDiscreteEventScheduler(TestCase):
         self.assertEqual(len(self.scheduler.event_queue), 1)
         event_time, _, event_agent = self.scheduler.event_queue[0]
         self.assertEqual(event_time, 5)
-        self.assertEqual(event_agent, self.agent1)
+        self.assertEqual(event_agent(), self.agent1)
 
     def test_schedule_event_with_float_time(self):
         self.scheduler.schedule_event(5.5, self.agent1)
         self.assertEqual(len(self.scheduler.event_queue), 1)
         event_time, _, event_agent = self.scheduler.event_queue[0]
         self.assertEqual(event_time, 5.5)
-        self.assertEqual(event_agent, self.agent1)
+        self.assertEqual(event_agent(), self.agent1)
 
     def test_schedule_in(self):
         self.scheduler.schedule_in(3, self.agent2)
         _, _, event_agent = self.scheduler.event_queue[0]
-        self.assertEqual(event_agent, self.agent2)
+        self.assertEqual(event_agent(), self.agent2)
         self.assertEqual(self.scheduler.get_next_event_time(), self.scheduler.time + 3)
 
     # Testing Event Execution and Time Advancement
@@ -337,6 +379,11 @@ class TestDiscreteEventScheduler(TestCase):
         with self.assertRaises(ValueError):
             self.scheduler.schedule_event(-1, self.agent1)
 
+    def test_invalid_aget_time(self):
+        with self.assertRaises(ValueError):
+            agent3 = MockAgent(3, self.model)
+            self.scheduler.schedule_event(2, agent3)
+
     def test_immediate_event_execution(self):
         # Current time of the scheduler
         current_time = self.scheduler.time
@@ -366,7 +413,7 @@ class TestDiscreteEventScheduler(TestCase):
         self.assertEqual(len(self.scheduler.event_queue), 1)
         event_time, _, event_agent = self.scheduler.event_queue[0]
         self.assertEqual(event_time, self.scheduler.time)
-        self.assertEqual(event_agent, new_agent)
+        self.assertEqual(event_agent(), new_agent)
 
         # Step the scheduler and check if the agent's step method is executed
         self.scheduler.step()