Generative Spaces (ABM)
In this workshop we will learn how to construct a ABM (Agent Based Model) with spatial behaviours, that is capable of configuring the space. This file is a simplified version of Generative Spatial Agent Based Models. For further information, you can find more advanced versions here:
PesudoCode
# initialization
import pandas, numpy, topogenesis, ...
define stencil & stencil_sq (neighborhood calculation)
# load all results
read lattice, preference, & adjacency matrix
order based on Area
adjust preferences to balance power effect
read dynamical results
fields = {}
for value in dunamical results:
fields[name] = value
# initialization agents
for a_id, a_prefs in preferences.iterrows():
update avail_index
a_eval = np.ones(len(avail_index))
# search for the best location inside available lattice
for f in fields:
a_eval *= fields[f][avail_index] ** a_prefs[f]
place agent to maximum point
update lattices
# calculate distance (estimate)
def distance(a_id,fns):
dist = []
for cen in fns_cens:
dist.append(distance between cen & the central point of a_id)
return dist
# the agent based model
while t < n_frames:
for a_id, a_prefs in preferences.iterrows():
# calculate free neighbors and neighbors for squareness
fns, fns_sq = [], []
neighs, neighs_sq = findneigh(stencil,loc[a_id]), findneigh(stencil_sq,loc[a_id])
for n in neighs:
if n is avail: fns.append(n)
for n in neighs_sq:
if n is avail: fns_sq.append(n)
# the evaluation process
if len(fns)>0:
a_eval = np.ones(len(fns))
# evaluate preference
for f in fields:
a_eval *= fields[f][fns] ** a_prefs[f]
# evaluate closeness to other agents
for s in agents:
a_eval *= distance(s,fns) ** weight
# evaluate squareness
fns_count = []
for fn in fns:
fns_count.append(fns_sq.count(fn))
a_eval *= fns_count ** square_weight
# if the agent has reached its desired space
calculate current_length
if current_length >= max_space:
(calculate preference & closeness)
# evaluate squareness
i_neighs_count = np.zeros(current_length)
for id,loc in enumerate(a_locs):
neighs = findneigh(stencil_sq,loc)
for n in neighs:
i_neighs_count[id] += (occ_lattice==a_id)[n]
i_eval *= i_neighs_count ** square_weight
# find lowest inner value
selected_int_inner = np.argmin(i_eval)
# find largerst fns value
selected_int = np.argmax(a_eval)
# if the agent has reached its desired space and the fns is better
if current_length >= max_space and min(a_eval) > max(i_eval):
remove the inner voxel
# add the new voxel
if not current_length >= max_space:
add selected new voxel
construct new_occ_lattice
frames.append(new_occ_lattice)
t += 1
# visualization and saving
visualize
save the frames to csv
0. Initialization
0.1. Load required libraries
# !pip install pyvista==0.28.1 ipyvtklink
import os
import topogenesis as tg
import pyvista as pv
import trimesh as tm
import pandas as pd
import numpy as np
import networkx as nx
from sklearn.cluster import KMeans
import scipy as sp
import pickle
import matplotlib.pyplot as plt
np.random.seed(0)
Parameter settings
# The number of frames of the growth model
n_frames = 1200
0.2. Define the Neighborhood (Stencil)
# creating neighborhood definition
stencil = tg.create_stencil("von_neumann", 1, 1)
# setting the center to zero
stencil.set_index([0,0,0], 0)
# creating neighborhood definition
stencil_sq = tg.create_stencil("von_neumann", 1, 1)
# setting the center to zero
# stencil_sq.set_index([0,0,0], 0)
stencil_sq.set_index([0,0,1], 0)
stencil_sq.set_index([0,0,-1], 0)
p = pv.Plotter(notebook=True)
grid = pv.UniformGrid()
grid.dimensions = np.array(stencil_sq.shape) + 1
grid.origin = [0,0,0]
grid.spacing = [1,1,1]
grid.cell_arrays["values"] = stencil_sq.flatten(order="F")
threshed = grid.threshold([0.9, 1.1])
p.add_mesh(threshed, show_edges=True, color="#ff8fa3", opacity=0.3)
p.show()
1. Setup the Environment
1.1. Load the envelope lattice as the avialbility lattice
# loading the lattice from csv
lattice_path = os.path.relpath('../Data/dynamic output/voxelized_envelope_cut.csv')
avail_lattice = tg.lattice_from_csv(lattice_path)
init_avail_lattice = tg.to_lattice(np.copy(avail_lattice), avail_lattice)
1.2 Load Program
# Load agent sizes
sizes_complete = pd.read_csv("../Data/raw data/Agent_sizes.csv")
Area = sizes_complete['Area']
sizes_complete = sizes_complete.sort_values(by = 'Area', ascending = 0)
display(sizes_complete)
# Calculating agent sizes in voxels
agent_areas = []
for i, row in sizes_complete.iterrows():
part = row['Area']
agent_areas.append(part)
# Load preference
square_weight = 0.1
program_complete = pd.read_csv("../Data/raw data/Programme_pref.csv")
program_complete = program_complete.drop(["sunlight"], 1)
program_complete['noise_field'] *= -1
program_complete['dist_entrance'] *= -0.5
program_complete['dist_fac'] *= -0.5
program_complete['Area'] = Area
program_complete = program_complete.sort_values(by = 'Area', ascending = 0)
program_complete = program_complete.drop(['Area'],1)
program_complete
# Simplify it
program_prefs = program_complete.drop(["space_name","space_id"], 1)
program_test = pd.read_csv("../Data/raw data/adjacency_matrix.csv")
program_test = program_test.drop(["Unnamed: 0"], 1)
1.2 Load the value fields
# loading the lattice from csv
fields = {}
for f in program_prefs.columns:
lattice_path = os.path.relpath('../Data/dynamic output/' + f + '.csv') # should care naming of csv
fields[f] = tg.lattice_from_csv(lattice_path)
1.3. Initialize the Agents
# initialize the occupation lattice
occ_lattice = avail_lattice * 0 - 1 # -1 means no agent is using it
agn_num = len(program_complete)
# Calculate best starting point
# The order to calculate is based on their sizes
agn_locs = [None]*agn_num
for a_id, a_prefs in program_complete.iterrows():
avail_index = np.array(np.where(avail_lattice)).T
a_eval = np.ones(len(avail_index))
for f in program_prefs.columns:
vals = fields[f][avail_index[:,0], avail_index[:,1], avail_index[:,2]]
# raise the the raw value to the power of preference weight of the agent
a_weighted_vals = vals ** a_prefs[f]
# multiply them to the previous weighted values
a_eval *= a_weighted_vals
for i in range(len(a_eval)):
if a_eval[i] == np.inf:
a_eval[i] = -1
selected_int = np.argmax(a_eval)
selected_ind = avail_index[selected_int]
n = 1
while True:
fns = avail_lattice.find_neighbours_masked(stencil, loc = selected_ind)
blocked = 0
for n in fns:
neigh_3d_id = np.unravel_index(n, avail_lattice.shape)
if occ_lattice[neigh_3d_id] != -1:
blocked += 1
if blocked >= 2:
selected_int = np.argsort(-a_eval,axis=0)[n]
selected_ind = avail_index[selected_int]
n += 1
else:
break
agn_locs[a_id]=[selected_ind]
avail_lattice[tuple(selected_ind)] = 0
occ_lattice[tuple(selected_ind)] = a_id
1.4. Visualize the environment
p = pv.Plotter(notebook=True)
# Set the grid dimensions: shape + 1 because we want to inject our values on the CELL data
grid = pv.UniformGrid()
grid.dimensions = np.array(occ_lattice.shape) + 1
# The bottom left corner of the data set
grid.origin = occ_lattice.minbound - occ_lattice.unit * 0.5
# These are the cell sizes along each axis
grid.spacing = occ_lattice.unit
# adding the boundingbox wireframe
p.add_mesh(grid.outline(), color="grey", label="Domain")
# adding axes
p.add_axes()
p.show_bounds(grid="back", location="back", color="#777777")
# making space list with index for the sargs
space_list = program_complete.get('space_name')
# formatting for the sarg annotation
space_list = space_list.to_dict()
sargs = dict(
shadow = True,
n_labels = 0,
italic = False,
fmt=" %.0f",
font_family="arial",
height = 0.6,
vertical = True,
position_x = 0.8,
position_y = 0.75)
# Add the data values to the cell data
grid.cell_arrays["Agents"] = occ_lattice.flatten(order="F").astype(int) # Flatten the array!
# filtering the voxels
threshed = grid.threshold([-0.1, agn_num - 0.9])
# adding the voxels
p.add_mesh(threshed, show_edges=True, opacity=1.0, show_scalar_bar=True, annotations = space_list, scalar_bar_args=sargs, cmap= "tab20")
# adding the availability lattice
# init_avail_lattice.fast_vis(p)
p.show()
2. ABM Simulation (Agent Based Space Occupation)
2.1. Function(s) for later use
# Calculate distance
# This is a stupid but doable estimate for faster calculation
lattice_cens = init_avail_lattice.centroids_threshold(-1)
def distance(a_id,fns):
fns_cens = []
for fn in fns:
fns_cens.append(lattice_cens[np.ravel_multi_index(fn,avail_lattice.shape)])
dist_m = []
for voxel_cen in fns_cens:
diff = voxel_cen - np.average(agn_locs[a_id],axis=0)
diff_p2 = diff**2
diff_p2s = diff_p2.sum()
dist = diff_p2s**0.5
dist_m.append(dist)
dist_m = np.array(dist_m)
return(dist_m)
2.2. Running the simulation
# make a deep copy of occupation lattice
cur_occ_lattice = tg.to_lattice(np.copy(occ_lattice), occ_lattice)
# initialzing the list of frames
frames = [cur_occ_lattice]
# setting the time variable to 0
t = 0
# main feedback loop of the simulation (for each time step ...)
while t<n_frames: # For every time
n_fns = 0
# for each agent ...
for a_id, a_prefs in program_complete.iterrows():
# retrieve the list of the locations of the current agent
a_locs = agn_locs[a_id]
# initialize the list of free neighbours
free_neighs = []
free_neighs_sq = []
# for each location of the agent
for loc in a_locs:
# retrieve the list of neighbours of the agent based on the stencil
neighs = avail_lattice.find_neighbours_masked(stencil, loc = loc)
neighs_sq = avail_lattice.find_neighbours_masked(stencil_sq, loc = loc)
# for each neighbour ...
for n in neighs:
# compute 3D index of neighbour
neigh_3d_id = np.unravel_index(n, avail_lattice.shape) # for checking availability
# if the neighbour is available...
if avail_lattice[neigh_3d_id]: # True -> Free
# add the neighbour to the list of free neighbours
free_neighs.append(neigh_3d_id) # Do not care about counting for the second time.
for n in neighs_sq:
neigh_3d_id = np.unravel_index(n, avail_lattice.shape)
if avail_lattice[neigh_3d_id]:
free_neighs_sq.append(neigh_3d_id)
# Investigate no possible free space
if not len(free_neighs)>0:
n_fns += 1
# check if found any free neighbour
if len(free_neighs)>0:
# convert free neighbours to a numpy array
fns = np.array(free_neighs)
# find the value of neighbours
# init the agent value array
a_eval = np.ones(len(fns)) # First set all values to 1
# This is the part to make amandment on choosing which block
# for each field...
for f in program_prefs.columns: # f represents names
# find the raw value of free neighbours...
vals = fields[f][fns[:,0], fns[:,1], fns[:,2]] # vals should be an array of length of fns
# raise the the raw value to the power of preference weight of the agent
a_weighted_vals = vals ** a_prefs[f]
# multiply them to the previous weighted values
a_eval *= a_weighted_vals
# now this part add the desirable space vicinity into the account
for s in program_test.columns:
s = int(s)
vals = distance(s,fns)
a_weighted_vals = vals ** program_test[str(a_id)][s]
a_eval *= a_weighted_vals
# This is the part that takes square shape of room into the account
free_neighs_count = []
# count can find how many times a location appears
for free_neigh in free_neighs:
free_neighs_count.append(free_neighs_sq.count(free_neigh))
# If a location is pointed for several times, it refers to more squared shape
a_weighted_square = np.array(free_neighs_count) ** square_weight
a_eval *= a_weighted_square
# Here, when it reaches the maximum space needed, start to evaluate inner voxels
current_length = np.copy(len(a_locs))
# This part means nothing, just for prevention of bugs
i_eval = np.zeros(current_length)
# max_space can be modified, or iterable (max_space[a_id])
max_space_raw = agent_areas[a_id] / (avail_lattice.unit[0] * avail_lattice.unit[1])
max_space = np.rint(max_space_raw)
# If the space occupied by this agent reaches designated
if current_length >= max_space:
i_eval = np.ones(current_length)
# The following parts do exactly the same calculations
# Calculating fields
for f in program_prefs.columns:
vals = fields[f][np.array(a_locs)[:,0], np.array(a_locs)[:,1], np.array(a_locs)[:,2]]
a_weighted_vals = vals ** a_prefs[f]
i_eval *= a_weighted_vals
# Calculating closeness
for s in program_test.columns:
s = int(s)
vals = distance(s,np.array(a_locs))
a_weighted_vals = vals ** program_test[str(a_id)][s]
i_eval *= a_weighted_vals
# Calculating squareness is bit different but the idea is same.
i_neighs_count = np.zeros(current_length)
for id,loc in enumerate(a_locs):
# Find neighbors of each point in location
neighs = init_avail_lattice.find_neighbours_masked(stencil_sq, loc = loc)
for n in neighs:
# Use occ_lattice == a_id to check if the neighborhoods of that point are the same agent.
neigh_3d_id = np.unravel_index(n, avail_lattice.shape)
i_neighs_count[id] += (occ_lattice==a_id)[neigh_3d_id]
i_weighted_square = np.array(i_neighs_count) ** square_weight
i_eval *= i_weighted_square
# select the inner with lowest evaluation
selected_int_inner = np.argmin(i_eval)
# select the neighbour with highest evaluation
selected_int = np.argmax(a_eval)
# If the agent reaches the required space, and the new free neighbour is better
# Then we have to remove the old voxel
if (current_length >= max_space) and i_eval[selected_int_inner] < a_eval[selected_int]:
selected_inner_3d_id = tuple(a_locs[selected_int_inner])
selected_inner_loc = a_locs[selected_int_inner]
agn_locs[a_id].pop(selected_int_inner)
avail_lattice[selected_inner_3d_id] = 1
occ_lattice[selected_inner_3d_id] = -1
# If the agent does not reach the required space, or we already removed old voxel
# Then we need to add new voxel
if current_length < max_space:
# find 3D integer index of selected neighbour
selected_neigh_3d_id = free_neighs[selected_int]
# find the location of the newly selected neighbour
selected_neigh_loc = np.array(selected_neigh_3d_id).flatten()
# add the newly selected neighbour location to agent locations
agn_locs[a_id].append(selected_neigh_loc)
# set the newly selected neighbour as UNavailable (0) in the availability lattice
avail_lattice[selected_neigh_3d_id] = 0
# set the newly selected neighbour as OCCUPIED by current agent
# (-1 means not-occupied so a_id)
occ_lattice[selected_neigh_3d_id] = a_id
# If the agent reaches the required space and the old voxels are better, we do nothing
# constructing the new lattice
new_occ_lattice = tg.to_lattice(np.copy(occ_lattice), occ_lattice)
# adding the new lattice to the list of frames
frames.append(new_occ_lattice) # For drawing graphes
# adding one to the time counter
print(t,"/",n_fns, end=" ")
t += 1
2.3. Visualizing the simulation
This is only for the sake of adjusting the model, can find a more detailed version at s6.2 of the repository.
p = pv.Plotter(notebook=True)
base_lattice = frames[0]
# Set the grid dimensions: shape + 1 because we want to inject our values on the CELL data
grid = pv.UniformGrid()
grid.dimensions = np.array(base_lattice.shape) + 1
# The bottom left corner of the data set
grid.origin = base_lattice.minbound - base_lattice.unit * 0.5
# These are the cell sizes along each axis
grid.spacing = base_lattice.unit
# adding the boundingbox wireframe
p.add_mesh(grid.outline(), color="grey", label="Domain")
# adding the availability lattice
init_avail_lattice.fast_vis(p)
# adding axes
p.add_axes()
p.show_bounds(grid="back", location="back", color="#aaaaaa")
def create_mesh(value):
f = int(value)
lattice = frames[f]
# Add the data values to the cell data
grid.cell_arrays["Agents"] = lattice.flatten(order="F").astype(int) # Flatten the array!
# filtering the voxels
threshed = grid.threshold([-0.1, agn_num - 0.9])
# adding the voxels
p.add_mesh(threshed, name='sphere', show_edges=True, opacity=1.0, show_scalar_bar=False)
return
p.add_slider_widget(create_mesh, [0, n_frames], title='Time', value=0, event_type="always", style="classic")
p.show(use_ipyvtk=True)
agn = 1
p = pv.Plotter(notebook=True)
base_lattice = frames[0]
# Set the grid dimensions: shape + 1 because we want to inject our values on the CELL data
grid = pv.UniformGrid()
grid.dimensions = np.array(base_lattice.shape) + 1
# The bottom left corner of the data set
grid.origin = base_lattice.minbound - base_lattice.unit * 0.5
# These are the cell sizes along each axis
grid.spacing = base_lattice.unit
# adding the boundingbox wireframe
p.add_mesh(grid.outline(), color="grey", label="Domain")
# adding the availability lattice
init_avail_lattice.fast_vis(p)
# adding axes
p.add_axes()
p.show_bounds(grid="back", location="back", color="#aaaaaa")
def create_mesh(value):
f = int(value)
lattice = frames[f]
# Add the data values to the cell data
grid.cell_arrays["Agents"] = lattice.flatten(order="F").astype(int) # Flatten the array!
# filtering the voxels
threshed = grid.threshold([agn-0.1, agn+0.9])
# adding the voxels
p.add_mesh(threshed, name='sphere', show_edges=True, opacity=1.0, show_scalar_bar=False)
return
p.add_slider_widget(create_mesh, [0, n_frames], title='Time', value=0, event_type="always", style="classic")
p.show(use_ipyvtk=True)
Saving to csv
for i, lattice in enumerate(frames):
csv_path = os.path.relpath('../Data/dynamic output/abm_animation/abm_f_'+ f'{i:03}' + '.csv')
lattice.to_csv(csv_path)
Credits
__author__ = "Shervin Azadi "
__license__ = "MIT"
__version__ = "1.0"
__url__ = "https://github.com/shervinazadi/spatial_computing_workshops"
__summary__ = "Spatial Computing Design Studio Workshop on Agent Based Models for Generative Spaces"