Seed Allocation and agent growth
In this notebook, we will first place the agents in an optimal voxel. We will do this by importing all of the previously computed lattices. When the agents have taken their place, they will grow in a specific direction and size.
The inputs needed for this notebook are the previously computed lattices, the program csv and the voxelized final envelope.
0. Initialization
0.1. Load required libraries
import os
import matplotlib as mpl
import matplotlib.pyplot as plt
from matplotlib import cm
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 copy
import scipy as sp
np.random.seed(0)
import pickle
# 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)
# 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 - center with 0 neighbours in xy plane
s_z = tg.create_stencil("von_neumann", 0, 1)
# setting the center to zero
s_z.set_index([0, 0, 0], 0)
# setting z neighbours to 1
s_z.set_index([0, 0,-1], 1)
s_z.set_index([0, 0, 1], 1)
# creating neighborhood definition - center with 0 neighbours in z plane
s_xy = tg.create_stencil("von_neumann", 1, 1)
# setting the center to zero
s_xy.set_index([0, 0, 0], 0)
# setting z neighbours to 0
s_xy.set_index([0, 0, 1], 0)
s_xy.set_index([0, 0, -1], 0)
0.3. Load the envelope lattice as the avialbility lattice
# loading the lattice from csv
context_path = os.path.relpath('../data/immediate_context.obj')
context_mesh = tm.load(context_path)
lattice_path = os.path.relpath('../data/final_envelope_new.csv')
avail_lattice = lattice_from_csv(lattice_path)
new_avail_lattice = (avail_lattice + 0.499) > 0.5
init_avail_lattice = tg.to_lattice(np.copy(avail_lattice), avail_lattice)
new_init_avail_lattice = tg.to_lattice(np.copy(new_avail_lattice), new_avail_lattice)
# visualizing the voxelized envelope
# convert trimesh to pv_mesh
def tri_to_pv(tri_mesh):
faces = np.pad(tri_mesh.faces, ((0, 0),(1,0)), 'constant', constant_values=3)
pv_mesh = pv.PolyData(tri_mesh.vertices, faces)
return pv_mesh
# initiating the plotter
p = pv.Plotter(notebook=True)
# fast visualization of the lattice
new_avail_lattice.fast_vis(p)
# plotting
p.show(use_ipyvtk=True)
0.4. Load Agents Information
# loading program (agents information) from CSV
# The program hase some extra collums, created for future expansion of the criteria of allocation and growth.
# We haven't been able to finish the implementation of those criteria:
prgm_path = os.path.relpath('../data/program.csv')
# start reading information after the first column and after the first row
agn_info = np.genfromtxt(prgm_path, delimiter=',')[1:, 1:]
# extract agent ids
agn_ids = agn_info[:, 0]
# extract agent preference to expand in the z-direction
agn_expandz = agn_info[:, 38]
# extract maximum voxels of each agent agent. This represents the maximum area & volume
agn_vox_req = agn_info[:, 39]
# extract criteria for silent level and noise repel
agn_silent_level = agn_info[:, 43]
agn_noise_repel = agn_info[:, 44]
# read the program with the use of panda's to create a dict. This is usefull to extract the right criteria for the right lattices.
agn_info_df = pd.read_csv('../data/program.csv')
agn_data = pd.read_csv(prgm_path)
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)
ent1_acc_path = os.path.relpath('../data/ent1_access.csv')
ent1_acc_lattice = lattice_from_csv(ent1_acc_path)
ent2_acc_path = os.path.relpath('../data/ent2_access.csv')
ent2_acc_lattice = lattice_from_csv(ent2_acc_path)
ent3_1_acc_path = os.path.relpath('../data/ent3.1_access.csv')
ent3_1_acc_lattice = lattice_from_csv(ent3_1_acc_path)
ent3_2_acc_path = os.path.relpath('../data/ent3.2_access.csv')
ent3_2_acc_lattice = lattice_from_csv(ent3_2_acc_path)
ent3_3_acc_path = os.path.relpath('../data/ent3.3_access.csv')
ent3_3_acc_lattice = lattice_from_csv(ent3_3_acc_path)
ent4_acc_path = os.path.relpath('../data/ent4_access.csv')
ent4_acc_lattice = lattice_from_csv(ent4_acc_path)
ent5_acc_path = os.path.relpath('../data/ent5_access.csv')
ent5_acc_lattice = lattice_from_csv(ent5_acc_path)
ent6_acc_path = os.path.relpath('../data/ent6_access.csv')
ent6_acc_lattice = lattice_from_csv(ent6_acc_path)
ent7_acc_path = os.path.relpath('../data/ent7_access.csv')
ent7_acc_lattice = lattice_from_csv(ent7_acc_path)
noise_acc_path = os.path.relpath('../data/HeerBokelweg_noise.csv')
noise_acc_lattice = lattice_from_csv(noise_acc_path)[:41, :21, :11]
public_green_acc_path = os.path.relpath('../data/public_greenery_highres.csv')
public_green_acc_lattice = lattice_from_csv(public_green_acc_path)
# this would be a lattice for the green corridor on the ground floor
# intern_green_acc_path = os.path.relpath('../data/green_openings_3.6.csv')
# intern_green_acc_lattice = lattice_from_csv(intern_green_acc_path)
# this would be a lattice to have functions attract to a facade facing the green corridor
# intern_facade_acc_path = os.path.relpath('../data/inner_facade_acces.csv')
# intern_facade_acc_lattice = lattice_from_csv(intern_facade_acc_path)
env_info_dict = {
"ent1_acc": ent1_acc_lattice ,
"ent2_acc": ent2_acc_lattice ,
"ent3.1_acc": ent3_1_acc_lattice ,
"ent3.2_acc": ent3_2_acc_lattice ,
"ent3.3_acc": ent3_3_acc_lattice ,
"ent4_acc": ent4_acc_lattice ,
"ent5_acc": ent5_acc_lattice ,
"ent6_acc": ent6_acc_lattice ,
"ent7_acc": ent7_acc_lattice ,
"sun_acc": sun_acc_lattice ,
"noise_acc": noise_acc_lattice ,
"public_green_acc": public_green_acc_lattice ,
#"intern_green_acc": intern_green_acc_lattice ,
#"inner_facade_acc": intern_facade_acc_lattice
}
env_info_list = []
env_info_dict_copy = env_info_dict. copy()
env_info_list.append(env_info_dict_copy)
env_info_base = copy.deepcopy(env_info_list)
# visualizing the info lattices
p = pv.Plotter(notebook=True)
info_val_list = list(env_info_dict.values())
info_key_list = env_info_dict.keys()
for i, k in enumerate(info_key_list):
print(i, k)
base_lattice = info_val_list[0]
# Create the spatial reference
grid = pv.UniformGrid()
# Set the grid dimensions: shape because we want to inject our values
grid.dimensions = base_lattice.shape
# The bottom left corner of the data set
grid.origin = base_lattice.minbound
# These are the cell sizes along each axis
grid.spacing = base_lattice.unit
p.add_mesh(tri_to_pv(context_mesh), style='surface')
def create_mesh(value):
f = int(value)
lattice = info_val_list[f]
# Add the data values to the cell data
grid.point_arrays["info"] = lattice.flatten(order="F") # Flatten the Lattice
# adding the volume
opacity = np.array([0,0.6,0.6,0.6,0.6,0.6,0.6])*1.5
p.add_volume(grid, cmap="coolwarm", name='sphere', clim=[0.0, 1.0],opacity=opacity, shade=True)
return
p.add_slider_widget(create_mesh, [0, len(info_val_list)-1], title='Time', value=0, event_type="always", style="classic", pointa=(0.1, 0.1), pointb=(0.9, 0.1))
p.show(use_ipyvtk=True)
# find the number of all voxels
vox_count = new_avail_lattice.size
# initialize the adjacency matrix
adj_mtrx = np.zeros((vox_count,vox_count))
# Finding the index of the available voxels in avail_lattice
avail_index = np.array(np.where(new_avail_lattice == 1)).T
# fill the adjacency matrix using the list of all neighbours
for vox_loc in avail_index:
# find the 1D id
vox_id = np.ravel_multi_index(vox_loc, new_avail_lattice.shape)
# retrieve the list of neighbours of the voxel based on the stencil
vox_neighs = new_avail_lattice.find_neighbours_masked(stencil, loc = vox_loc)
# iterating over the neighbours
for neigh in vox_neighs:
# setting the entry to one
adj_mtrx[vox_id, neigh] = 1.0
# construct the graph
g = nx.from_numpy_array(adj_mtrx)
# compute the distance of all voxels to all voxels using floyd warshal algorithm#
# this will be done once to save computation time, therefore the next line is commented. When the envelope changes, the caluclation has to be done again
# dist_mtrx = nx.floyd_warshall_numpy(g)
# this will save the calculated dist_mtrx. Again this will be done once since we calculate the dist_mtrx once.
# pickle.dump(dist_mtrx, open( "../data/dist_mtrx_final.p", "wb" ) )
# the saved dist_mtrx is loaded in a variable to use
dist_mtrx_loaded = pickle.load( open( "../data/dist_mtrx_final.p", "rb" ) )
1. ABM Simulation
1.1. Initialize the Agents
# initialize the occupation lattice
occ_lattice = new_avail_lattice * 0 - 1
# Finding the index of the available voxels in avail_lattice
avail_flat = new_avail_lattice.flatten()
avail_index = np.array(np.where(new_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 = []
# retrieving the entrance access value of the free neighbours
for a_id in agn_ids:
voxel_vals = []
pot_voxels = []
# retrieve agent preferences;
a_pref_dict = agn_info_df.loc[a_id]
# Voxel Evaluation Loop
for pot_vox in avail_index:
if new_avail_lattice[tuple(pot_vox)]:
global_vox_value = 1.0
# for every lattice in the environment informations
for info_key, info_lattice in env_info_dict.items():
# Here we utilise Fuzzy Logics to be able to compare different layers
# of environmental information and evaluate the voxel for the agent.
vox_val = info_lattice[tuple(pot_vox)]
# agn_vox_val = np.power(vox_val, a_pref[i])
agn_vox_val = np.power(vox_val, a_pref_dict[info_key])
global_vox_value *= agn_vox_val
# add the neighbour value to the list of values
voxel_vals.append(global_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_id = tuple(pot_voxels[selected_int].T)
# 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.append([selected_neigh_loc])
# set the newly selected neighbour as UNavailable (0) in the availability lattice
new_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
def dynamic_noise_lattice(agn_locs, new_avail_lattice):
# define the noise range
noise_range = [10.0, 60.0]
# initialize noise sources
noise_src_points = []
noise_src_levels = []
# iterate over agents
for a_id in range(len(agn_locs)):
# extract agent locations
a_locs = agn_locs[a_id]
# retrieve the silent level of the agent
a_noise_level_mapped = 1 - agn_data.loc[a_id]["silent_level"]
# mapping the [0,1] values to noise level (db)
a_noise_level = a_noise_level_mapped * (noise_range[1] - noise_range[0]) + noise_range[0]
# for each agent location
for a_loc in a_locs:
# append the noise source information
noise_src_points.append(a_loc)
noise_src_levels.append(a_noise_level)
# convert to numpy array
noise_src_points = np.array(noise_src_points)
# create full lattice
full_lattice = new_avail_lattice * 0 + 1
# extract the coordiantes of the centroid of all voxel
vox_centroids = full_lattice.centroids
# extract voxel indices of all voxels
vox_indices = np.array(np.where(full_lattice==1)).T
# initializing the sum lattice of noise
sum_noise_lats = new_avail_lattice * 0.0
# for each source of noise
for src_point, src_level in zip(noise_src_points,noise_src_levels):
# initialize the occupation lattice
dist_latice = new_avail_lattice * 0.0
for cen, ind in zip(vox_centroids, vox_indices):
# compute the euclidian distance
dist_latice[tuple(ind)] = sp.spatial.distance.euclidean(cen, src_point)
# computing the noise lattice from dist lattice
noise_latice = src_level - 20 * np.log10(dist_latice) - 8
# summing
sum_noise_lats += np.power(10, noise_latice / 10.0)
# computing the final aggregation
agg_noise_lats = 10 * np.log10(sum_noise_lats)
# normalizing the noise values
normalized_silence_lattice = 1 - (agg_noise_lats - np.min(agg_noise_lats)) / (np.max(agg_noise_lats) - np.min(agg_noise_lats))
return normalized_silence_lattice
1.2. Running the simulation
dynamic_info = {"noise_repel": dynamic_noise_lattice}
# 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
# choosing the number of frames. Our maximum of voxels is 620, therefore if you want to see the full growth choose 620.
# 620 frames takes a while to calculate, choose something around 10 to quick test the code
n_frames = 620
# Simulation Loop
# main feedback loop of the simulation (for each time step ...)
while t<n_frames:
# this lets functions repel each other. This does take very long to calculate.
# for info_key, info_function in dynamic_info.items():
# env_info_dict[info_key] = info_function(agn_locs, avail_lattice)
# this part is about the dynamic relations
agn_acc_list = []
# for each agent
for a_id in range(agn_num):
# find which voxel occupied by destination function (roof garden)
a_locs = agn_locs[a_id]
# organize 3d indices vertically
a_locs_array = np.array(a_locs).T
# extract 1d index for all the locations of this agent
a_locs_1d = np.ravel_multi_index(a_locs_array, new_avail_lattice.shape)
# extract distance field to those voxels from the distance matrix
voxel_dists = dist_mtrx_loaded[a_locs_1d]
# combine dist fields of all voxels of this agent
combined_dists = voxel_dists.min(axis=0)
max_valid = np.ma.masked_invalid(combined_dists).max()
# set the infinities to one more than the maximum valid values
combined_dists[combined_dists == np.inf] = max_valid + 1
# mapping the values from (0, max) to (1, 0)
mapped_dists = 1 - combined_dists / max_valid
# constructing the lattice
agn_acc_lattice = tg.to_lattice(mapped_dists.reshape(new_avail_lattice.shape), new_avail_lattice)
agn_acc_list.append(agn_acc_lattice)
# combine agn_dist fields with base env info
env_info = env_info_base + agn_acc_list
# Agent Loop
# for each agent ...
for a_id in range(agn_num):
# retrieve the list of the locations of the current agent
a_locs = agn_locs[a_id]
# retrieve the list of the prefernece to grow in the z plane of the current agent
a_expandz = agn_expandz[int(a_id)]
# retrieve the list of the maximum amount of voxels for each function
a_vox_req = agn_vox_req[a_id]
# retrieve the list of the amount of voxels of each function
a_vox_count = len(a_locs)
# initialize the list of free neighbours
free_neighs = []
free_neighs_1d = []
free_neigh_weights = []
# Location loop
# for each location of the agent
for loc in a_locs:
# here is decided if a function grows in the xy plane or in the xyz plane
# retrieve the list of neighbours of the agent based on the stencil
z_neighs = new_avail_lattice.find_neighbours_masked(s_z, loc = loc)
# retrieve the list of neighbours of the agent based on the stencil
xy_neighs = new_avail_lattice.find_neighbours_masked(s_xy, loc = loc)
# neighs are stacked in sequence horizontally
neighs = np.hstack([z_neighs,xy_neighs])
# neigh_weights are stacked in sequence horizontally. z preference is variable, xy prefernce is 1.0.
# because it is 1.0 is will always grow horizontally and it will also grow vertically based on the weight of a_expandz
neigh_weights = np.hstack([np.full(z_neighs.shape, a_expandz), np.full(xy_neighs.shape,1.0)])
# for each neighbour ...
# compare neighs and neigh_weights
for n, nw in zip(neighs, neigh_weights):
# compute 3D index of neighbour
neigh_3d_id = np.unravel_index(n, new_avail_lattice.shape)
# if the neighbour is available...
if new_avail_lattice[neigh_3d_id]:
# add the neighbour to the list of free neighbours
free_neighs.append(neigh_3d_id)
free_neigh_weights.append(nw)
free_neighs_1d.append(n)
# check if found any free neighbour
if len(free_neighs)>0:
# convert free neighbours to a numpy array
free_neighs = np.array(free_neighs)
free_neighs_1d = np.array(free_neighs_1d)
# uniq_neigh: neighbours that occur more than once are only called once
# unique_to_norm_order: numbers in free_neighs_1d replaced by positions of unig_neigh
# unique_neighs_count: counts the amount of times a neighbour appears in free_neighs_1d
uniq_neighs, unique_to_norm_order, unique_neighs_count = np.unique(free_neighs_1d, return_counts=True, return_inverse=True)
# voxel value is raised when a voxel appears more often in free_neighs_1d
# 1.5 is the indicator for preference to grow more in a square shape choosing a reoccuring neighbour.
# this could be a variable number for each space
extra_weight = (unique_neighs_count - 1) * 1.5 + 1
# retrieving the entrance access value of the free neighbours
neigh_vals = []
# retrieve agent preferences
a_pref_dict = agn_info_df.loc[a_id]
# Neighbour Evaluation Loop
for neigh, neigh_w in zip(free_neighs,free_neigh_weights):
neigh_value = neigh_w
# for every lattice in the environment informations
# for i, info_lattice in enumerate(env_info):
for info_key, info_lattice in env_info_dict.items():
# Here we utilise Fuzzy Logics to be able to compare different layers
# of environmental information and evaluate the voxel for the agent.
vox_val = info_lattice[tuple(neigh)]
agn_vox_val = np.power(vox_val, a_pref_dict[info_key])
neigh_value *= agn_vox_val
# add the neighbour value to the list of values
neigh_vals.append(neigh_value)
# convert to numpy array
neigh_vals = np.array(neigh_vals)
# extra weight is added to voxel values
neigh_vals *= extra_weight[unique_to_norm_order]
# select the neighbour with highest value
selected_int = np.argmax(neigh_vals)
selected_neigh_val = neigh_vals[selected_int]
# find 3D intiger index of selected neighbour
selected_neigh_3d_id = tuple(free_neighs[selected_int].T)
# find the location of the newly selected neighbour
selected_neigh_loc = np.array(selected_neigh_3d_id).flatten()
# here we check if a function can still grow based on its current amount of voxels and the maximum required
if a_vox_count < a_vox_req:
# 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
new_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 required amount of voxels has been reached, no more new neighbours will be added to the location list of that agent
else:
pass
# 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
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 avilability lattice
# new_init_avail_lattice.fast_vis(p)
# new_avail_lattice.fast_vis(p)
p.add_mesh(tri_to_pv(context_mesh), style='surface')
# Make a dictionary for the annotations
annotations = {
0: "ent1",
1: "ent2",
2: "ent3.1",
3: "ent3.2",
4: "ent3.3",
5: "ent4",
6: "ent5",
7: "ent6",
8: "ent7",
9: "student",
10: "assisted",
11: "starter",
12: "restaurant",
13: "shop",
14: "cocooking",
15: "pub",
16: "gym",
17: "arcade",
18: "cinema",
19: "office",
20: "cowork",
21: "library",
22: "fablabs",
23: "catering",
24: "catering2",
25: "catering3",
26: "coffeehub"
}
# make a dictionary for customizations
sargs = dict(
title_font_size=14,
label_font_size=6,
shadow=True,
n_labels=0,
italic=False,
fmt="%.0f",
font_family="arial",
height=0.8,
vertical=True,
position_x=1.0,
position_y=0.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=annotations, scalar_bar_args=sargs, cmap="tab20")
return
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
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)
Credits
__author__ = "Shervin Azadi and Pirouz Nourian
__changes_made_by__ "Xam Adan"
__license__ = "MIT"
__version__ = "1.0"
__url__ = "https://github.com/frankvahstal/spatial_computing_workshops"
__summary__ = "Spatial Computing Design Studio Workshop on MCDA and Path Finding for Generative Spatial Relations"