Generative Relations: MCDA script
Enabling agents to utilize MCDA (Multi Criteria Decision Analyses) in their spatial behaviors.
input:
Lattice info; CSV files from the sun/sky/quitness access and entrance and/or specific parts of the lattice accesses
CSV with the table of values for each lattice info (as described before) and remaining necessary values such as (maz z coordinate, max voxel count, stencil id and evaluation per agent) included within the space names and space id's.
output:
Lattice frames (growth) saved as CSV.
0. Initialization
0.1. Load required libraries
import os
import topogenesis as tg
import pyvista as pv
import trimesh as tm
import numpy as np
import networkx as nx
import pandas as pd
import scipy as sp
import pickle
import matplotlib.pyplot as plt
np.random.seed(0)
# extra import function
def lattice_from_csv(file_path):
# read metadata
meta_df = pd.read_csv(file_path, nrows=3)
shape = np.array(meta_df['shape'])
unit = np.array(meta_df['unit'])
minbound = np.array(meta_df['minbound'])
# read lattice
lattice_df = pd.read_csv(file_path, skiprows=5)
# create the buffer
buffer = np.array(lattice_df['value']).reshape(shape)
# create the lattice
l = tg.to_lattice(buffer, minbound=minbound, unit=unit)
return l
0.2. Define the Neighborhood (Stencil)
0.2.1. Basic stencils with z axes differences
# creating neighborhood definition for stencil that is 1.8m high
s_1 = tg.create_stencil("von_neumann", 1, 1)
# setting the center to zero
s_1.set_index([0,0,0], 0)
# creating neighborhood definition for stencil that is 3.6m high
s_2 = tg.create_stencil("von_neumann", 1, 2)
# setting the center to zero
s_2.set_index([0, 0, 0], 0)
s_2.set_index([0, 0, 1], 0)
s_2.set_index([0, 0, 2], 1)
s_2.set_index([0, 0,-1], 0)
s_2.set_index([0, 0,-2], 1)
# setting the center to zero
s_2.set_index([0,0,0], 0)
# creating neighborhood definition for stencil that is 5.4m high
s_3 = tg.create_stencil("von_neumann", 1, 3)
# setting the center to zero
s_3.set_index([0, 0, 0], 0)
s_3.set_index([0, 0, 1], 0)
s_3.set_index([0, 0, 3], 1)
s_3.set_index([0, 0,-1], 0)
s_3.set_index([0, 0,-3], 1)
# setting the center to zero
s_3.set_index([0,0,0], 0)
# creating neighborhood definition for stencil that is 7.2m high
s_4 = tg.create_stencil("von_neumann", 1, 4)
# setting the center to zero
s_4.set_index([0, 0, 0], 0)
s_4.set_index([0, 0, 1], 0)
s_4.set_index([0, 0, 4], 1)
s_4.set_index([0, 0,-1], 0)
s_4.set_index([0, 0,-4], 1)
# setting the center to zero
s_4.set_index([0,0,0], 0)
# creating neighborhood definition for stencil that is 9m high
s_5 = tg.create_stencil("von_neumann", 1, 5)
# setting the center to zero
s_5.set_index([0, 0, 0], 0)
s_5.set_index([0, 0, 1], 0)
s_5.set_index([0, 0, 5], 1)
s_5.set_index([0, 0,-1], 0)
s_5.set_index([0, 0,-5], 1)
# setting the center to zero
s_5.set_index([0,0,0], 0)
# listing the stencils in order to make them correspond later with the spaces and their height requirement
stencils = [s_1, s_2, s_3, s_4, s_5]
0.2.2. Visualization stencils
0.2.3. Listing all stencils
stencil = s_5
# initiating the plotter
p = pv.Plotter(notebook=True)
# Create the spatial reference
grid = pv.UniformGrid()
# Set the grid dimensions: shape because we want to inject our values
grid.dimensions = np.array(stencil.shape) + 1
# The bottom left corner of the data set
grid.origin = [0,0,0]
# These are the cell sizes along each axis
grid.spacing = [5,5,5]
# Add the data values to the cell data
grid.cell_arrays["values"] = stencil.flatten(order="F") # Flatten the stencil
threshed = grid.threshold([0.9, 1.1])
# adding the voxels: light red
p.add_mesh(threshed, show_edges=True, color="white", opacity=0.3)
# plotting
p.show(use_ipyvtk=True)
0.2.4 Stencil visualization (for checking)
0.3. Load the envelope lattice as the avialbility lattice
# loading the lattice from csv
lattice_path = os.path.relpath('../data/highres_envelope.csv')
avail_lattice = lattice_from_csv(lattice_path)
init_avail_lattice = tg.to_lattice(np.copy(avail_lattice), avail_lattice)
# loading the lattice from csv
lattice_path = os.path.relpath('../data/avail_lattice_good_voxels.csv')
avail_lattice_good_voxels = lattice_from_csv(lattice_path)
init_avail_lattice_good_voxels = tg.to_lattice(np.copy(avail_lattice_good_voxels), avail_lattice_good_voxels)
avail_lattice*= avail_lattice_good_voxels
0.4. Load Agents Information
# loading program (agents information) from CSV
prgm_path = os.path.relpath('../data/program_exported.csv')
agn_info = pd.read_csv('../data/program_exported.csv',delimiter=";")
agn_ids = agn_info["space_id"].values
agn_prefs = agn_info
agn_prefs
0.5. Initialize environment information layers from Sun Access Lattice and Entrance Access Lattice
# loading the lattice from csv
sun_acc_path = os.path.relpath('../data/sun_access_highres.csv')
sun_acc_lattice = lattice_from_csv(sun_acc_path)
ent_acc_highres_public_path = os.path.relpath('../data/ent_access_highres_public.csv')
ent_acc_public_lattice = lattice_from_csv(ent_acc_highres_public_path)
ent_acc_highres_housing_path = os.path.relpath('../data/ent_access_highres_housing.csv')
ent_acc_housing_lattice = lattice_from_csv(ent_acc_highres_housing_path)
ent_acc_highres_gym_path = os.path.relpath('../data/ent_access_highres_gym.csv')
ent_acc_gym_lattice = lattice_from_csv(ent_acc_highres_gym_path)
ent_acc_highres_parking_path = os.path.relpath('../data/ent_access_highres_parking.csv')
ent_acc_parking_lattice = lattice_from_csv(ent_acc_highres_parking_path)
ent_acc_highres_comcen_path = os.path.relpath('../data/ent_access_highres_comcen.csv')
ent_acc_comcen_lattice = lattice_from_csv(ent_acc_highres_comcen_path)
highres_sky_acc_path = os.path.relpath('../data/sky_access_highres.csv')
sky_acc_lattice = lattice_from_csv(highres_sky_acc_path)
highres_quietness_acc_path = os.path.relpath('../data/quietness_highres.csv')
quietness_acc_lattice = lattice_from_csv(highres_quietness_acc_path)
groundfloor_acc_path = os.path.relpath('../data/ent_access_highres_groundfloor.csv')
groundfloor_acc_lattice = lattice_from_csv(groundfloor_acc_path)
# list the environment information layers (lattices)
env_info = {"sun_acc": sun_acc_lattice + 0.001,
"ent_acc_public": ent_acc_public_lattice + 0.001,
"ent_acc_housing": ent_acc_housing_lattice + 0.001,
"ent_acc_gym": ent_acc_gym_lattice + 0.001,
"ent_acc_parking": ent_acc_parking_lattice + 0.001,
"ent_acc_comcen": ent_acc_comcen_lattice + 0.001,
"sky_acc": sky_acc_lattice + 0.001,
"quietness_acc": quietness_acc_lattice + 0.001,
"ground_floor_acc": groundfloor_acc_lattice + 0.001}
# defining other factors in csv
# defining stencil id
stencil_id = stencils
# area to use in simulation
room_area = []
# check if all the csv files loaded are highres
for key, info_lattice in env_info.items():
print(key, info_lattice.shape)
1. ABM Simulation
Basic run
Making defs
# Function for checking the availability (Since it is repeated several times in the main loop)
def check_avail(avail_lattice, ind, a_stencil_id):
condition = 1
ind_array = np.array(ind)
for step in range(a_stencil_id + 1):
new_ind_array = ind_array + np.array([0,0,step])
condition *= avail_lattice[tuple(new_ind_array)]
return condition
# Function for evaluation (Since it is repeated several times in the main loop)
def eval_voxel(vox, env_info, a_pref):
global_vox_value = 1.0
# for every lattice in the environment informations
for key, info_lattice in env_info.items():
# Here we utilise Fuzzy Logics to be able to compare different layers
# of environmental information and evaluate the voxel for the agent.
# This method is introduced, and generalised in Pirouz Nourian dissertation:
# section 5.7.3, pp. 201-208, eq. 57. You can refer to this section for
# comprehensive mathematical details.
vox_val = info_lattice[tuple(vox)]
agn_vox_val = np.power(vox_val, a_pref[key])
global_vox_value *= agn_vox_val
return global_vox_value
# Function for the occupation and departure (Since it is repeated several times in the main loop)
def mult_occupation(selected_neigh_3d_address, a_id, a_height, agn_locs, agn_src_locs, occ_lattice, avail_lattice, departure=False):
# Doing this for x times in the z axis with x coming from a_height
for step in range(a_height):
#giving a step to the regular occupation in order to run this for every step in the z aces
new_address = selected_neigh_3d_address + np.array([0,0,step])
# check if there's enough space in the z axis
if new_address[2] < occ_lattice.shape[2]:
# make tuple of the address
selected_neigh_3d_id = tuple(new_address)
# find the location of the newly selected neighbour
selected_neigh_loc = np.array(selected_neigh_3d_id).flatten()
if departure==False:
# add the newly selected neighbour location to agent locations
agn_locs[a_id].append(selected_neigh_loc)
if step == 0:
agn_src_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
else:
# remove the newly selected neighbour location to agent locations
a_locs_list = [list(loc) for loc in agn_locs[a_id]]
try:
ind = a_locs_list.index(list(selected_neigh_loc))
agn_locs[a_id].pop(ind)
# set the selected neighbour as available (1) in the availability lattice
avail_lattice[selected_neigh_3d_id] = 1
# set the newly selected neighbour as NOT OCCUPIED by current agent
occ_lattice[selected_neigh_3d_id] = -1
except:
pass
return (agn_locs, agn_src_locs, occ_lattice, avail_lattice)
# Simple version of finding neighbours provided by Shervin in order to fix the growth issue (agents didnt grow on the bottom level of the lattice).
# This function takes a stencil, a lattice and the address of a voxel in that lattice and it returns the index of the neighbours of that voxel
def find_neighbours_masked(lattice, stencil, loc):
neigh_locs = np.argwhere(stencil) - stencil.origin + loc
neigh_filter = np.all(neigh_locs > 0, axis=1) * np.all(neigh_locs < np.array(lattice.shape), axis=1)
return(neigh_locs[neigh_filter])
Intialization
# initialize the occupation lattice
occ_lattice = avail_lattice * 0 - 1
# Finding the index of the available voxels in avail_lattice
avail_flat = avail_lattice.flatten()
avail_index = np.array(np.where(avail_lattice == 1)).T
# count the number of spaces (rows) and intiialize an agent for each space
agn_num = len(agn_info)
# adding the origins to the agents locations
agn_locs = [[] for a_id in agn_ids]
agn_src_locs = [[] for a_id in agn_ids]
agn_upper = []
# retrieving the entrance access value of the free neighbours
for a_id in agn_ids:
voxel_vals = []
pot_voxels = []
# retrieve agent preferences
a_pref = agn_prefs.loc[a_id]
a_stencil_id = agn_prefs["stencil_id"][a_id]
stencil = stencils[a_stencil_id]
# max height (z dimension)
a_max_z = agn_prefs["max_z"][a_id]
# use avail_index
# Voxel Evaluation Loop
for pot_vox in avail_index:
if check_avail(avail_lattice, tuple(pot_vox), a_stencil_id) and pot_vox[2] < a_max_z:
# eval voxel
vox_value = eval_voxel(pot_vox, env_info, a_pref)
# add the neighbour value to the list of values
voxel_vals.append(vox_value)
pot_voxels.append(pot_vox)
# convert to numpy array
voxel_vals = np.array(voxel_vals)
# convert to numpy array
pot_voxels = np.array(pot_voxels)
# select the neighbour with highest value
selected_int = np.argmax(voxel_vals)
# find 3D intiger index of selected neighbour
selected_neigh_3d_address = tuple(pot_voxels[selected_int].T)
# find the location of the newly selected neighbour
agn_origins = np.array(selected_neigh_3d_address).flatten()
agn_locs, agn_src_locs, occ_lattice, avail_lattice = mult_occupation(selected_neigh_3d_address,
a_id,
a_stencil_id + 1,
agn_locs,
agn_src_locs,
occ_lattice,
avail_lattice)
1.2. Running the simulation
Basic run
# empty lists to fill in order to track satisfaction and area growth
satisfaction_record = []
area_record = []
# 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
n_frames = 2500
# check size of used csv (reach of the code)
print(agn_info.shape)
# max voxel count per space (here we do the stencil type + 1
#times the amount of area needed in order to obtain the total amount of voxels that need to be occupied by the script)
# pick room area data
a_room_vox = agn_prefs["room_area"]
# pick stencil id's and do +1 because we start with 0 instead of 1 (id 2 would have stencil 3 high otherwise)
a_room_stencil = agn_prefs["stencil_id"] + 1
# obtain the max amount of voxels needing to be occupied by doing the room area times the height of the space (stencil)
a_room_voxels = a_room_stencil * a_room_vox
# print in order to check below if you obtain correct values
print(a_room_voxels)
# Simulation Loop
# main feedback loop of the simulation (for each time step ...)
while t<n_frames:
all_agn_sat = []
agn_loc_vals = []
all_agn_area = []
# for each agent ... (evaluation/satisfaction)
for a_id in range(agn_num):
# creating variable that takes per id
a_locs = agn_locs[a_id]
a_src_locs = agn_src_locs[a_id]
# retrieving the entrance access value of the free neighbours
loc_vals = []
# retrieve agent preferences
a_pref = agn_prefs.loc[a_id]
# Neighbour Evaluation Loop
for loc in a_src_locs:
# eval voxel
loc_value = eval_voxel(tuple(loc), env_info, a_pref)
loc_vals.append(loc_value)
agn_loc_vals.append(loc_vals)
agent_satisfaction = np.mean(loc_value)
all_agn_sat.append(agent_satisfaction)
# end of the satisfaction for loop
satisfaction_record.append(all_agn_sat)
# Main Agent Loop
for a_id in range(agn_num):
# creating variable that takes per id
a_locs = agn_locs[a_id]
a_src_locs = agn_src_locs[a_id]
# creating variable that takes stencil per id
a_stencil_id = agn_prefs["stencil_id"][a_id]
# making a variable that gives the height of the stencils in voxel (coincidentally +1 since all stencils grow with 1 voxel per id)
a_height = a_stencil_id + 1
# variable for stencil values
stencil = stencils[a_stencil_id]
# making a variable of the result per id (for the loop)
a_room_voxel = a_room_voxels[a_id]
# initialize the list of free neighbours
free_neighs = []
# Location loop
# for each location of the agent
for loc in a_src_locs:
# retrieve the list of neighbours of the agent based on the stencil
neighs_3d = find_neighbours_masked(avail_lattice, s_2, loc=loc)
# for each neighbour ...
for neigh_3d_id in neighs_3d:
if check_avail(avail_lattice, neigh_3d_id, a_stencil_id) and neigh_3d_id[2] < a_max_z:
# add the neighbour to the list of free neighbours
free_neighs.append(neigh_3d_id)
# check if found any free neighbour
if len(free_neighs)>0:
# convert free neighbours to a numpy array
free_neighs = np.array(free_neighs)
# retrieving the entrance access value of the free neighbours
neigh_vals = []
# retrieve agent preferences
a_pref = agn_prefs.loc[a_id]
# Neighbour Evaluation Loop
for neigh in free_neighs:
# vox eval
neigh_value = eval_voxel(tuple(neigh), env_info, a_pref)
neigh_vals.append(neigh_value)
# convert to numpy array
neigh_vals = np.array(neigh_vals)
# select the neighbour with highest value
selected_neighbour = np.argmax(neigh_vals)
# find 3D intiger index of selected neighbour
selected_neigh_3d_address = free_neighs[selected_neighbour].T
agent_satisfaction = all_agn_sat[a_id]
# removing for better voxels
a_evaluation = a_pref["evaluation"]
# normal growth condition
if len(agn_locs[a_id]) < a_room_voxel:
# Ocupation
agn_locs, agn_src_locs, occ_lattice, avail_lattice = mult_occupation(selected_neigh_3d_address,
a_id,
a_height,
agn_locs,
agn_src_locs,
occ_lattice,
avail_lattice)
# if evaluation value is higher than the agent satisfaction
elif a_evaluation > agent_satisfaction:
# select the lowest valued voxel for that agent
loc_vals = agn_loc_vals[a_id]
selected_loc = np.argmin(loc_vals)
# Departure
if neigh_vals[selected_neighbour] > loc_vals[selected_loc]:
# Running stencil departure with the function
vox_loc = agn_src_locs[a_id][selected_loc]
agn_locs, agn_src_locs, occ_lattice, avail_lattice = mult_occupation(vox_loc,
a_id,
a_height,
agn_locs,
agn_src_locs,
occ_lattice,
avail_lattice,
departure=True)
# Ocupation after departure with the function
agn_locs, agn_src_locs, occ_lattice, avail_lattice = mult_occupation(free_neighs[selected_neighbour],
a_id,
a_height,
agn_locs,
agn_src_locs,
occ_lattice,
avail_lattice)
# end of swapping condition
# end of dynamic condition
# end of available free neighbour condition
# inside the agent iteration loop
all_agn_area.append(len(agn_locs[a_id]))
# end of the agent iteration loop
area_record.append(all_agn_area)
# 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)
# adding one to the time counter
t += 1
# satisfaction tracking
satisfaction_df = pd.DataFrame.from_records(satisfaction_record, columns=["Agent_" + str(a_id) for a_id in range(agn_num)])
satisfaction_df.index = ["Frame_" + str(f) for f in range(len(satisfaction_record))]
satisfaction_df
#save the table
table = satisfaction_df.to_csv('../data/2500frames.csv')
#load the table
table = pd.read_csv('../data/2500frames.csv')
# dividing the table
sat1 = table[["Agent_0", "Agent_1", "Agent_2", "Agent_3", "Agent_4", "Agent_5", "Agent_6",
"Agent_7", "Agent_8", "Agent_9", "Agent_10", "Agent_12", "Agent_13", "Agent_14", "Agent_15",
"Agent_18"]]
# set limit
sat1 = sat1[0:100]
# dividing the table
sat2 = table[["Agent_11", "Agent_17", "Agent_19"]]
# set limit
sat2 = sat2[0:800]
# dividing the table
sat3 = table[["Agent_16"]]
# set limit
sat3 = sat3[0:2500]
#plot graph
df = sat1
df.plot(legend=True, figsize=(15, 10), xlabel="new x", ylabel="new y", title = "xx", colormap="tab20");
plt.legend(bbox_to_anchor=(1.0, 0.5))
#plot graph
df = sat2
df.plot(legend=True, figsize=(15, 10), xlabel="new x", ylabel="new y", title = "xx", colormap="tab10");
plt.legend(bbox_to_anchor=(1.0, 0.5))
#plot graph
df = sat3
df.plot(legend=True, figsize=(15, 10), xlabel="new x", ylabel="new y", title = "xx", colormap="tab10");
plt.legend(bbox_to_anchor=(1.0, 0.5))
#area tracking
spacesize_df = pd.DataFrame.from_records(area_record, columns=["Agent_" + str(a_id) for a_id in range(agn_num)])
spacesize_df.index = ["Frame_" + str(f) for f in range(len(area_record))]
spacesize_df
#saving the tables to compare if script is runned again with different results
pickle.dump(frames, open("../data/frames.p", "wb" ))
#saving the tables to compare if script is runned again with different results
frames = pickle.load(open("../data/frames.p", "rb" ))
for frame in frames:
frame.bounds = occ_lattice.bounds
frame.unit = occ_lattice.unit
1.3. Visualizing the simulation
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 avilability lattice
init_avail_lattice_good_voxels.fast_vis(p)
# adding axes
p.add_axes()
p.show_bounds(grid="back", location="back", color="#aaaaaa")
# making space list with index for the sargs
space_list = agn_prefs.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 = 1.05,
position_y = 1)
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=True, annotations = space_list, scalar_bar_args=sargs, cmap="nipy_spectral")
return
# new cell in order to maintain vis size
p.add_slider_widget(create_mesh, [0, n_frames], title='Time', value=0, event_type="always", style="classic", pointa=(0.1, 0.1), pointb=(0.9, 0.1))
p.show(use_ipyvtk=True)
2.3. Saving lattice frames in CSV
# save for every frame
for i, lattice in enumerate(frames):
csv_path = os.path.relpath('../data/abm_mcda/abm_f_'+ f'{i:03}' + '.csv')
lattice.to_csv(csv_path)
#or just end frame
csv_path = os.path.relpath('../data/final_lattice.csv')
lattice.to_csv(csv_path)
# An issue occured regarding specific agents not occuring when loading in and visualizing the lattice,
#it was solved this way:
#load in final lattice
lattice_path = os.path.relpath('../data/final_lattice.csv')
avail_lattice_good_voxels = lattice_from_csv(lattice_path)
init_avail_lattice_good_voxels = tg.to_lattice(np.copy(avail_lattice_good_voxels), avail_lattice_good_voxels)
#instead of values, make it true/false
high_threshold = 50
low_threshold = -0.9
new_avail_lattice = ((init_avail_lattice_good_voxels < high_threshold) * (init_avail_lattice_good_voxels > low_threshold))
#save it
csv_path = os.path.relpath('../data/fixed_final_lattice.csv')
new_avail_lattice.to_csv(csv_path)
Credits
__author__ = "Shervin Azadi and Pirouz Nourian"
_Editor_ = "Eda Akaltun"
__license__ = "MIT"
__version__ = "1.0"
__url__ = "https://github.com/shervinazadi/spatial_computing_workshops"
__summary__ = "Spatial Computing Design Studio Workshop on MCDA and Path Finding for Generative Spatial Relations"