This is the repost of the following question as suggested by @HoboProber .
Again, if you don't know what is Schelling's model of segregation, you can read it here.
The Schelling model of segregation is an agent-based model that illustrates how individual tendencies regarding neighbors can lead to segregation. In the Schelling model, agents occupy cells of rectangular space. A cell can be occupied by a single agent only. Agents belong to one of two groups and are able to relocate according to the fraction of friends (i.e., agents of their own group) within a neighborhood around their location. The model's basic assumption is as follows: an agent, located in the center of a neighborhood where the fraction of friends f is less than a predefined tolerance threshold F (i.e., f < F), will try to relocate to a neighborhood for which the fraction of friends is at least f (i.e., f ≥ F)
You can get the shapefile here if you want to try it.
I have implemented most of the changes suggested and also did split up my code in 4 different classes. Here is my complete code.
geo_schelling_populate.py
import numpy as np
from shapely.geometry import Point
import geopandas as gp
import pandas as pd
class geo_schelling_populate:
""" Generate the coordinates in a polygon (In this case a map of the state)
on the basis of the given spacing and then randomly assign coordiantes
to different races and as empty houses. It also takes ratio and
empty_ratio in consideration.
Parameters
----------
shapefile : str
It is a string pointing to a geospatial vector data file which
has information for an american state's geometry data.
spacing : float
It defines the distance between the points in coordinate units.
empty_ratio : float
What percent of houses need to be kept empty.
ratio : float
What is the real ratio between the majority and minority in
reality in a given state.
ratio : int
Number of races. Currently the model is tested on 2 races, but
maybe in future support fot more races will be added.
Attributes
----------
shapefile : str
Same as in parameter section
spacing : float
Same as in parameter section
empty_ratio : float
Same as in parameter section
demographic_ratio : float
Same as ratio in the parameter section
races : int
Same as in parameter section
shape_file : geopandas.GeoDataFrame
Pandas DataFrame with a geometry column generated from the .shp
file
df_allhouses : pandas.DataFrame
Pandas DataFrame which contains all the different coordinates
generated inside a specific state.
df_emptyhouses : pandas.DataFrame
Pandas DataFrame which contains all the coordiantes which are
emptyhouses.
df_agenthouses : pandas.DataFrame
Pandas DataFrame which contains all the coordinates associated
with a race.
Methods
-------
_generate_points()
Private function! It generates the coordinates inside a given
shape. Returns a dataframe with all the coordiantes generated.
populate()
It populates the coordinates with either a race or denotes it
as an empty house. Returns a tuple with pandas dataframe of
empty houses and agent houses. Again agents are races in this
case.
"""
def __init__(
self,
shapefile,
spacing,
empty_ratio,
ratio,
races=2,
):
self.shapefile = shapefile
self.spacing = spacing
self.empty_ratio = empty_ratio
self.demographic_ratio = ratio
self.races = races
self.shape_file = \
gp.read_file(r"%s"%shapefile).explode().reset_index().drop(columns=['level_0'
, 'level_1'])
self.df_allhouses = pd.DataFrame([])
self.df_emptyhouses = pd.DataFrame([])
self.df_agenthouses = pd.DataFrame([])
def _generate_points(self, polygon):
"""It returns a DataFrame with all the coordiantes inside a certain
shape passed in as an parameter.
Parameters
----------
polygon : shapely.geometry.Polygon
A polygon object which contains the geometry of a county in
a state.
Returns
-------
A pandas DataFrame with all the coordiantes generated inside the
polygon object.
"""
(minx, miny, maxx, maxy) = polygon.bounds
x_coords = np.arange(np.floor(minx), int(np.ceil(maxx)),
self.spacing)
y_coords = np.arange(np.floor(miny), int(np.ceil(maxy)),
self.spacing)
grid = np.column_stack((np.meshgrid(x_coords,
y_coords)[0].flatten(),
np.meshgrid(x_coords,
y_coords)[1].flatten()))
df_points = pd.DataFrame.from_records(grid, columns=['X', 'Y'])
df_points = df_points[df_points[['X', 'Y']].apply(lambda x: \
Point(x[0], x[1]).within(polygon),
axis=1)]
return df_points.round(2)
def populate(self):
""" Populates the coordinates by assigning them a certain race or by
assigning them as empty houses.
Parameters
----------
No parameters
Returns
-------
A tuple which consist two pandas DataFrames. One contains
the coordiantes who have empty houses and the other contains
the coordiante with a race assigned to them.
"""
pd.set_option('mode.chained_assignment', None)
self.df_allhouses = \
pd.concat(iter(self.shape_file.geometry.apply(lambda x: \
self._generate_points(x))), ignore_index=True)
empty_ratio_seperator = round(self.empty_ratio
* len(self.df_allhouses))
self.df_allhouses['Agent/Empty'] = \
np.random.permutation(np.asarray(['Empty']
* empty_ratio_seperator
+ ['Agent'] * (len(self.df_allhouses)
- empty_ratio_seperator)))
self.df_emptyhouses = \
self.df_allhouses[self.df_allhouses['Agent/Empty']
== 'Empty'][['X', 'Y'
]].reset_index(drop=True)
self.df_agenthouses = \
self.df_allhouses[self.df_allhouses['Agent/Empty']
== 'Agent'][['X', 'Y'
]].reset_index(drop=True)
demographic_ratio_seperator = round(self.demographic_ratio
* len(self.df_agenthouses))
self.df_agenthouses['Race'] = \
np.random.permutation(np.asarray([1]
* demographic_ratio_seperator
+ [2] * (len(self.df_agenthouses)
- demographic_ratio_seperator)))
return (self.df_emptyhouses, self.df_agenthouses)
geo_schelling_update.py
import numpy as np
import math
class geo_schelling_update:
"""Updates the position of the races and the empty houses on the basis of
similarity threshold. Agents change positions if they have enough agents
from their own races in a certain neighbourhood. In this scenario,
neighbourhood is defined as the eight nearest coordinate around the
coordinate we are checking satisfaction for.
Parameters
----------
n_iterations : int
Maximum number of times, update method should run.
spacing : float
It defines the distance between the points in coordinate units.
np_agenthouses : numpy array
df_agenthouses converted to the numpy array for faster
computations.
np_emptyhouses : numpy array
df_emptyhouses converted to the numpy array for faster
computations.
similarity_threshold : float
What percent of similarity people want from their neighbours.
Attributes
----------
spacing : float
Same as in parameter section
n_iterations : int
Same as in parameter section
np_agenthouses : numpy array
Same as in parameter section
np_emptyhouses : numpy array
Same as in parameter section
similarity_threshold : float
Same as in parameter section
Methods
-------
_is_unsatisfied()
Private function! It checks if an agent is satisfied or not at
a certain position.
_move_to_empty()
Private function! Moves an unsatisified agent to a random empty
house.
_update_helper()
Private functions! A helper function to help update method with
numpy's apply_along_axis method.
update()
It updates array as long as it reaches maximum number of
iterations or reaches a point where no agent is unsatisfied.
Getter
------
get_agenthouses()
Returns an array with all the agents at a certain coordinate.
"""
def __init__(
self,
n_iterations,
spacing,
np_agenthouses,
np_emptyhouses,
similarity_threshold,
):
self.spacing = spacing
self.n_iterations = n_iterations
self.np_emptyhouses = np_emptyhouses
self.np_agenthouses = np_agenthouses
self.similarity_threshold = similarity_threshold
def _is_unsatisfied(self, x, y):
""" Checks if an agent is unsatisfied at a certain position.
Parameters
----------
x : float
x coordinate of the agent being checked
y : float
y coordinate of the agent being checked
Returns
-------
True or False based on if the agent is satisfied or not.
"""
race = np.extract(np.logical_and(np.equal(self.np_agenthouses[:
, 0], x), np.equal(self.np_agenthouses[:, 1],
y)), self.np_agenthouses[:, 2])[0]
euclid_distance1 = round(math.hypot(self.spacing,
self.spacing), 4)
euclid_distance2 = self.spacing
total_agents = \
np.extract(np.logical_or(np.equal(np.round(np.hypot(self.np_agenthouses[:
, 0] - x, self.np_agenthouses[:, 1] - y), 4),
euclid_distance1),
np.equal(np.round(np.hypot(self.np_agenthouses[:
, 0] - x, self.np_agenthouses[:, 1] - y), 4),
euclid_distance2)), self.np_agenthouses[:, 2])
if total_agents.size == 0:
return False
else:
return total_agents[total_agents == race].size \
/ total_agents.size < self.similarity_threshold
def _move_to_empty(self, x, y):
"""Moves the agent to a new position if the agent is unsatisfied.
Parameters
----------
x : float
x coordinate of the agent being checked
y : float
y coordinate of the agent being checked
Returns
-------
None
"""
race = np.extract(np.logical_and(np.equal(self.np_agenthouses[:
, 0], x), np.equal(self.np_agenthouses[:, 1],
y)), self.np_agenthouses[:, 2])[0]
(x_new, y_new) = \
self.np_emptyhouses[np.random.choice(self.np_emptyhouses.shape[0],
1), :][0]
self.np_agenthouses = \
self.np_agenthouses[~np.logical_and(self.np_agenthouses[:,
0] == x, self.np_agenthouses[:, 1]
== y)]
self.np_agenthouses = np.vstack([self.np_agenthouses, [x_new,
y_new, race]])
self.np_emptyhouses = \
self.np_emptyhouses[~np.logical_and(self.np_emptyhouses[:,
0] == x_new, self.np_emptyhouses[:, 1]
== y_new)]
self.np_emptyhouses = np.vstack([self.np_emptyhouses, [x, y]])
def _update_helper(self, agent):
"""Helps the update function with number of changes made in every
iterations.
Parameters
----------
agent : tuple
x and y coordinates for the agent's position.
Returns
-------
1 if the position of the agent is changed, 0 if not.
"""
if self._is_unsatisfied(agent[0], agent[1]):
self._move_to_empty(agent[0], agent[1])
return 1
else:
return 0
def update(self):
"""Main player in updating the array with all the agents' position.
It updates the array until it reaches the iteration limit or the
number of changes become 0.
Parameters
----------
None
Returns
-------
None
"""
for i in np.arange(self.n_iterations):
np_oldagenthouses = self.np_agenthouses.copy()
n_changes = np.sum(np.apply_along_axis(self._update_helper,
1, np_oldagenthouses))
print('n Changes---->' + str(n_changes))
print(i)
if n_changes == 0:
break
def get_agenthouses(self):
return self.np_agenthouses
geo_schelling_data.py
from shapely.geometry import Point
import pandas as pd
class geo_schelling_data:
"""Get the important data from the simulation required for the further
analysis. Following data is obtained:
1) County Name
2) Majority Population
3) Minority Population
4) Total Population
5) Percentage of Majority Population
6) Percentage of Minority Population
Parameters
----------
np_agenthouses : numpy array
df_agenthouses converted to the numpy array for faster
computations.
pd_shapefile : Pandas DataFrame
Pandas DataFrame with a geometry column generated from the .shp
file.
Attributes
----------
np_agenthouses : numpy array
Same as in parameter section
pd_shapefile : Pandas DataFrame
Same as in parameter section
df_county_data : Pandas DataFrame
Pandas DataFrame with important data mentioned above.
Methods
-------
_get_number_by_county()
Private Function! It returns the number of majority and minority
in a polygon.
get_county_data()
Returns the pandas DataFrame with the important data required for
the further analysis as mentioned above.
"""
def __init__(self, np_agenthouses, pd_shapefile):
self.np_agenthouses = np_agenthouses
self.pd_shapefile = pd_shapefile
self.df_county_data = pd.DataFrame([])
def _get_number_by_county(self, geometry):
"""Returns the number of minority agents and majority agents within a
polygon.
Parameters
----------
geometry : shapely.geometry.Polygon
A polygon object which contains the geometry of a county in
a state.
Returns
-------
A tuple with of minority agents and majority agents
"""
num_majority = len([x[2] for x in list(self.np_agenthouses)
if Point(x[0], x[1]).within(geometry)
and x[2] == 1.0])
num_minority = len([x[2] for x in list(self.np_agenthouses)
if Point(x[0], x[1]).within(geometry)
and x[2] == 2.0])
return (num_majority, num_minority)
def get_county_data(self):
"""Does calculation on Minority and Majority agents' number and returns
a pandas dataframe for further analysis.
Parameters
----------
None
Returns
-------
Pandas DataFrame
"""
self.df_county_data['CountyName'] = self.pd_shapefile.NAMELSAD
self.df_county_data['MajorityPop'] = \
self.pd_shapefile.geometry.apply(lambda x: \
self._get_number_by_county(x)[0])
self.df_county_data['MinorityPop'] = \
self.pd_shapefile.geometry.apply(lambda x: \
self._get_number_by_county(x)[1])
self.df_county_data['TotalPop'] = \
self.df_county_data['MajorityPop'] \
+ self.df_county_data['MinorityPop']
self.df_county_data['MajorityPopPercent'] = \
self.df_county_data[['TotalPop', 'MajorityPop'
]].apply(lambda x: (0 if x['TotalPop']
== 0 else x['MajorityPop'] / x['TotalPop']), axis=1)
self.df_county_data['MinorityPopPercent'] = \
self.df_county_data[['TotalPop', 'MinorityPop'
]].apply(lambda x: (0 if x['TotalPop']
== 0 else x['MinorityPop'] / x['TotalPop']), axis=1)
return self.df_county_data
geo_schelling_plot.py
from matplotlib import pyplot as plt
class geo_schelling_plot:
"""Visualize the simulation as it helps us interpret and analyze the
results.
Parameters
----------
np_agenthouses : numpy array
df_agenthouses converted to the numpy array for faster
computations.
pd_shapefile : Pandas Series
Pandas Series of the geometry column generated from the .shp
file.
Attributes
----------
np_agenthouses : numpy array
Same as in parameter section
pd_shapefile : Pandas Series
Same as in parameter section
Methods
-------
plot()
Plots the agents and shape on the same graph.
"""
def __init__(self, np_agenthouses, pd_geometry):
self.np_agenthouses = np_agenthouses
self.pd_geometry = pd_geometry
def plot(self):
"""Plots the agents and shape on the same graph so that one can see
where the different agents lies. It is very helpful in the case when
the simulation is finished and you can see the changes occured in
the plot.
Yellow is majority
Violet is minority
Parameters
----------
None
Returns
-------
None
"""
fig, ax = plt.subplots(figsize=(10, 10))
self.pd_geometry.plot(ax=ax, color='white', edgecolor='black',linewidth=4)
ax.scatter(self.np_agenthouses[:, 0], self.np_agenthouses[:,
1], c=self.np_agenthouses[:, 2])
ax.set_xticks([])
ax.set_yticks([])
run.py
from geo_optimized.geo_schelling_data import geo_schelling_data
from geo_optimized.geo_schelling_plot import geo_schelling_plot
from geo_optimized.geo_schelling_populate import geo_schelling_populate
from geo_optimized.geo_schelling_update import geo_schelling_update
shapefile = "Path to the shapefile"
spacing = 0.05
empty_ratio = 0.2
similarity_threshhold = 0.65
n_iterations = 100
ratio = 0.41
if __name__ == '__main__':
schelling_populate = geo_schelling_populate(shapefile,
spacing,
empty_ratio,
ratio)
df_emptyhouses, df_agenthouses = schelling_populate.populate()
np_agenthouses = df_agenthouses.to_numpy()
np_emptyhouses = df_emptyhouses.to_numpy()
pd_geometry = schelling_populate.shape_file.geometry
schelling_update = geo_schelling_update(n_iterations,
spacing,
np_agenthouses,
np_emptyhouses,similarity_threshhold)
schelling_update.update()
np_agenthouses = schelling_update.get_agenthouses()
pd_shapefile = schelling_populate.shape_file
schelling_plot = geo_schelling_plot(np_agenthouses,pd_geometry)
schelling_plot.plot()
schelling_get_data = geo_schelling_data(np_agenthouses,pd_shapefile)
df_county_data = schelling_get_data.get_county_data()
I got rid of most of the for loops and replace the code with pandas and numpy to do some faster looping and also some vectorize operations.
I believe this code is much more cleaner and faster than the previous one but still some improvements in speed and documentation can be made.
I believe I can still make some operations vectorize but don't know how to do it. If someone can suggest those that would be great. Also if someone can help me with the better documentation of the code, that would be great too. Currently I am using numpy style docstrings.
Also if someone can help me optimize my code with numba to achieve C level speed that would be great. I believe that geo_schelling_update.py can be speed up with numba but I am unable to do so.
.shp. There are some associated files which are needed gis.stackexchange.com/a/262509 \$\endgroup\$ – Maarten Fabré Sep 18 at 12:39