diff --git a/agent.py b/agent.py index c0e7b1d..48e6e7c 100644 --- a/agent.py +++ b/agent.py @@ -5,6 +5,12 @@ This model implements the actual agents on the grid (a.k.a. the ants) License: AGPL 3 (see end of file) (C) Alexander Bocken, Viviane Fahrni, Grace Kagho +""" + +""" +TO DISCUSS: +Is the separation of energy and sensitivity useful? + """ import numpy as np import numpy.typing as npt @@ -12,14 +18,10 @@ from mesa.agent import Agent from mesa.space import Coordinate class RandomWalkerAnt(Agent): - def __init__(self, unique_id, model, look_for_pheromone=None, - energy_0=1, - pheromone_drop_rate_0 : dict[str, float]={"A": 80, "B": 80}, - sensitivity_0=0.99, - alpha=0.6, drop_pheromone=None, - betas : dict[str, float]={"A": 0.0512, "B": 0.0512}, - sensitivity_decay_rate=0.01, - sensitivity_max = 300 + def __init__(self, unique_id, model, + look_for_pheromone=None, + drop_pheromone=None, + sensitivity_max = 30000, ) -> None: super().__init__(unique_id=unique_id, model=model) @@ -27,18 +29,12 @@ class RandomWalkerAnt(Agent): self._next_pos : None | Coordinate = None self._prev_pos : None | Coordinate = None - self.look_for_pheromone = look_for_pheromone - self.drop_pheromone = drop_pheromone - self.energy = energy_0 #TODO: use - self.sensitivity_0 = sensitivity_0 - self.sensitivity = self.sensitivity_0 - self.pheromone_drop_rate = pheromone_drop_rate_0 - self.alpha = alpha + self.look_for_pheromone : str|None = look_for_pheromone + self.drop_pheromone : str|None = drop_pheromone + self.energy : float = self.model.e_0 + self.sensitivity : float = self.model.s_0 + self.pheromone_drop_rate : float = self.model.q_0 self.sensitivity_max = sensitivity_max - self.sensitivity_decay_rate = sensitivity_decay_rate - self.betas = betas - self.threshold : dict[str, float] = {"A": 0, "B": 0} - def sens_adj(self, props, key) -> npt.NDArray[np.float_] | float: """ @@ -68,7 +64,7 @@ class RandomWalkerAnt(Agent): # TODO: proper nonlinear response, not just clamping if props > self.sensitivity_max: return self.sensitivity_max - if props > self.threshold[key]: + if props > self.model.q_tr: return props else: return 0 @@ -78,9 +74,29 @@ class RandomWalkerAnt(Agent): arr.append(self.sens_adj(prop, key)) return np.array(arr) + def _get_resistance_weights(self, positions=None): + if positions is None: + positions = self.neighbors() + # bit round-about but self.model.grid.fields['res'][positions] + # gets interpreted as slices, not multiple singular positions + resistance = np.array([ self.model.grid.fields['res'][x,y] for x,y in positions ]) + easiness = np.max(self.model.grid.fields['res']) - resistance + 1e-15 # + epsilon to not divide by zero + weights = easiness/ np.sum(easiness) + + return weights + def _choose_next_pos(self): + def _pick_from_remaining_five(remaining_five): + """ + """ + weights = self._get_resistance_weights(remaining_five) + random_index = np.random.choice(range(len(remaining_five)), p=weights) + self._next_pos = remaining_five[random_index] + self._prev_pos = self.pos + if self._prev_pos is None: - i = np.random.choice(range(6)) + weights = self._get_resistance_weights() + i = np.random.choice(range(6),p=weights) assert(self.pos is not self.neighbors()[i]) self._next_pos = self.neighbors()[i] self._prev_pos = self.pos @@ -89,71 +105,94 @@ class RandomWalkerAnt(Agent): if self.searching_food: for neighbor in self.front_neighbors: if self.model.grid.is_food(neighbor): + self.model.grid.fields['food'][neighbor] -= 1 # eat + #resets + self.pheromone_drop_rate = self.model.q_0 + self.sensitivity = self.model.s_0 + self.energy = self.model.e_0 + + #now look for other pheromone self.drop_pheromone = "B" self.look_for_pheromone = "A" - self.sensitivity = self.sensitivity_0 self._prev_pos = neighbor self._next_pos = self.pos + return elif self.searching_nest: for neighbor in self.front_neighbors: if self.model.grid.is_nest(neighbor): + #resets + self.pheromone_drop_rate = self.model.q_0 + self.sensitivity = self.model.s_0 + self.energy = self.model.e_0 + self.look_for_pheromone = "A" # Is this a correct interpretation? self.drop_pheromone = "A" - self.sensitivity = self.sensitivity_0 self._prev_pos = neighbor self._next_pos = self.pos # recruit new ants - for agent_id in self.model.get_unique_ids(self.model.num_new_recruits): - if self.model.schedule.get_agent_count() < self.model.num_max_agents: + for agent_id in self.model.get_unique_ids(self.model.N_r): + if self.model.schedule.get_agent_count() < self.model.N_m: agent = RandomWalkerAnt(unique_id=agent_id, model=self.model, look_for_pheromone="B", drop_pheromone="A") agent._next_pos = self.pos self.model.schedule.add(agent) self.model.grid.place_agent(agent, pos=neighbor) + return - # follow positive gradient + # follow positive gradient with likelihood self.sensitivity if self.look_for_pheromone is not None: + # Calculate gradient front_concentration = [self.model.grid.fields[self.look_for_pheromone][cell] for cell in self.front_neighbors ] front_concentration = self.sens_adj(front_concentration, self.look_for_pheromone) current_pos_concentration = self.sens_adj(self.model.grid.fields[self.look_for_pheromone][self.pos], self.look_for_pheromone) gradient = front_concentration - np.repeat(current_pos_concentration, 3).astype(np.float_) - # TODO: if two or more neighbors have same concentration randomize? Should be unlikely with floats though + index = np.argmax(gradient) if gradient[index] > 0: - self._next_pos = self.front_neighbors[index] - self._prev_pos = self.pos + # follow positive gradient with likelihood self.sensitivity + p = np.random.uniform() + if p < self.sensitivity: + self._next_pos = self.front_neighbors[index] + self._prev_pos = self.pos + else: + other_neighbors = self.neighbors().copy() + other_neighbors.remove(self.front_neighbors[index]) + _pick_from_remaining_five(other_neighbors) return # do biased random walk p = np.random.uniform() - if p < self.alpha: + # TODO: This completely neglects resistance, relevant? + if p < self.model.alpha: self._next_pos = self.front_neighbor self._prev_pos = self.pos else: - # need copy() as we would otherwise remove the tuple from all possible lists (aka python "magic") other_neighbors = self.neighbors().copy() other_neighbors.remove(self.front_neighbor) - random_index = np.random.choice(range(len(other_neighbors))) - self._next_pos = other_neighbors[random_index] - self._prev_pos = self.pos + _pick_from_remaining_five(other_neighbors) def step(self): - self.sensitivity -= self.sensitivity_decay_rate - self._choose_next_pos() - self._adjust_pheromone_drop_rate() + self.sensitivity -= self.model.d_s + self.energy -= self.model.grid.fields['res'][self.pos] * self.model.d_e + # Die and get removed if no energy + if self.energy < self.model.e_min: + self.model.schedule.remove(self) + else: + self._choose_next_pos() + self._adjust_pheromone_drop_rate() def _adjust_pheromone_drop_rate(self): if(self.drop_pheromone is not None): - self.pheromone_drop_rate[self.drop_pheromone] -= self.pheromone_drop_rate[self.drop_pheromone] * self.betas[self.drop_pheromone] + self.pheromone_drop_rate -= self.pheromone_drop_rate * self.model.beta def drop_pheromones(self) -> None: # should only be called in advance() as we do not use hidden fields if self.drop_pheromone is not None: - self.model.grid.fields[self.drop_pheromone][self.pos] += self.pheromone_drop_rate[self.drop_pheromone] + self.model.grid.fields[self.drop_pheromone][self.pos] += self.pheromone_drop_rate def advance(self) -> None: self.drop_pheromones() diff --git a/hexplot.py b/hexplot.py new file mode 100755 index 0000000..62eb3ca --- /dev/null +++ b/hexplot.py @@ -0,0 +1,31 @@ +#!/bin/python +import numpy as np +import matplotlib.pyplot as plt + +def plot_hexagon(A, title=None): + X, Y = np.meshgrid(range(A.shape[0]), range(A.shape[-1])) + X, Y = X*2, Y*2 + + # Turn this into a hexagonal grid + for i, k in enumerate(X): + if i % 2 == 1: + X[i] += 1 + Y[:,i] += 1 + + fig, ax = plt.subplots() + im = ax.hexbin( + X.reshape(-1), + Y.reshape(-1), + C=A.reshape(-1), + gridsize=int(A.shape[0]/2) + ) + + # the rest of the code is adjustable for best output + ax.set_aspect(1) + ax.set(xlim=(-4, X.max()+4,), ylim=(-4, Y.max()+4)) + ax.axis(False) + plt.colorbar(im) + + if(title is not None): + plt.title(title) + plt.show(block=False) diff --git a/main.py b/main.py index c9c3c15..0d780ce 100755 --- a/main.py +++ b/main.py @@ -6,6 +6,7 @@ execute via `python main.py` in terminal or only UNIX: `./main.py` License: AGPL 3 (see end of file) (C) Alexander Bocken, Viviane Fahrni, Grace Kagho """ +import array from model import ActiveWalkerModel from agent import RandomWalkerAnt import numpy as np @@ -16,10 +17,11 @@ from mesa.datacollection import DataCollector from multihex import MultiHexGrid def main(): - check_pheromone_exponential_decay() - check_ant_sensitivity_linear_decay() - check_ant_pheromone_exponential_decay() - check_ants_follow_gradient() + pass + # check_pheromone_exponential_decay() + # check_ant_sensitivity_linear_decay() + # check_ant_pheromone_exponential_decay() + # check_ants_follow_gradient() def check_pheromone_exponential_decay(): """ @@ -107,6 +109,7 @@ def check_ant_pheromone_exponential_decay(): num_initial_roamers = 1 num_max_agents = 100 nest_position : Coordinate = (width //2, height //2) + num_food_sources = 0; max_steps = 1000 model = ActiveWalkerModel(width=width, height=height, @@ -179,8 +182,27 @@ def check_ants_follow_gradient(): model.step() -if __name__ == "__main__": - main() +# if __name__ == "__main__": +# main() + +from model import kwargs_paper_setup1 as kwargs +model = ActiveWalkerModel(**kwargs) + +from hexplot import plot_hexagon +# a = np.zeros_like(model.grid.fields['food']) +# a[np.nonzero(model.grid.fields['food'])] = 1 +# plot_hexagon(a, title="Nest locations") +# plot_hexagon(model.grid.fields['res'], title="Resistance Map") + + +from tqdm import tqdm as progress_bar +for _ in progress_bar(range(model.max_steps)): + model.step() +# agent_densities = model.datacollector.get_model_vars_dataframe()["agent_dens"] +# mean_dens = np.mean(agent_densities) +# norm_dens = mean_dens/np.max(mean_dens) +# plot_hexagon(norm_dens, title="Ant density overall") +# plt.show() diff --git a/model.py b/model.py index 10cee4a..ae29610 100644 --- a/model.py +++ b/model.py @@ -16,40 +16,131 @@ from mesa.time import SimultaneousActivation from mesa.datacollection import DataCollector from agent import RandomWalkerAnt +kwargs_paper_setup1 = { + "width": 100, + "height": 100, + "N_0": 20, + "N_m": 100, + "N_r": 5, + "alpha": 0.6, + "gamma": 0.001, + "beta": 0.0512, + "d_s": 0.001, + "d_e": 0.001, + "s_0": 0.99, + "e_0": 0.99, + "q_0": 80, + "q_tr": 1, + "e_min": 0, + "nest_position": (49,49), + "N_f": 5, + "food_size" : 55, + "max_steps": 8000, + "resistance_map_type" : None, +} + +kwargs_paper_setup2 = { + "width": 100, + "height": 100, + "N_0": 20, + "N_m": 100, + "N_r": 5, + "alpha": 0.6, + "gamma": 0.01, + "beta": 0.0512, + "d_s": 0.001, + "d_e": 0.001, + "s_0": 0.99, + "e_0": 0.99, + "q_0": 80, + "q_tr": 1, + "e_min": 0, + "nest_position": (49,49), + "N_f": 5, + "food_size" : 550, + "max_steps": 8000, + "resistance_map_type" : None, +} + class ActiveWalkerModel(Model): - def __init__(self, width : int, height : int , num_max_agents : int, - num_initial_roamers : int, + def __init__(self, width : int, height : int, + N_0 : int, # number of initial roamers + N_m : int, # max number of ants + N_r : int, # number of new recruits + alpha : float, #biased random walk + beta : float, # decay rate drop rate + gamma : float, # decay rate pheromone concentration fields + d_s : float, # decay rate sensitvity + d_e : float, # decay rate energy + s_0 : float, # sensitvity reset + e_0 : float, # energy reset + q_0 : float, # initial pheromone level + q_tr : float, # threshold under which ant cannot distinguish concentrations + e_min : float, # energy at which walker dies nest_position : Coordinate, - num_food_sources=5, - food_size=10, + N_f=5, #num food sources + food_size= 55, max_steps:int=1000, + resistance_map_type=None, ) -> None: super().__init__() - fields=["A", "B", "nests", "food"] + + self.N_m : int = N_m # max number of ants + self.N_r : int = N_r # number of new recruits + self.alpha : float = alpha # biased random walk if no gradient + self.gamma : float = gamma # decay rate pheromone concentration fields + self.beta : float = beta # decay rate drop rate + self.d_s : float = d_s # decay rate sensitvity + self.d_e : float = d_e # decay rate energy (get's multiplied with resistance) + self.s_0 : float = s_0 # sensitvity reset + self.e_0 : float = e_0 # energy reset + self.q_0 : float = q_0 # pheromone drop rate reset + self.q_tr : float = q_tr # threshold under which ant cannot distinguish concentrations + self.e_min : float = e_min # energy at which walker dies + self.N_f : int = N_f #num food sources + + fields=["A", "B", "nests", "food", "res"] self.schedule = SimultaneousActivation(self) self.grid = MultiHexGridScalarFields(width=width, height=height, torus=True, fields=fields) + + if resistance_map_type is None: + self.grid.fields["res"] = np.ones((width, height)).astype(float) + elif resistance_map_type == "perlin": + # perlin generates isotropic noise which may or may not be a good choice + # pip3 install git+https://github.com/pvigier/perlin-numpy + from perlin_numpy import ( + generate_fractal_noise_2d, + generate_perlin_noise_2d, + ) + noise = generate_perlin_noise_2d(shape=(width,height), res=((10,10))) + normalized_noise = (noise - np.min(noise))/(np.max(noise) - np.min(noise)) + self.grid.fields["res"] = normalized_noise + else: + # possible other noise types: simplex or value + raise NotImplemented(f"{resistance_map_type=} is not implemented.") + + + + self._unique_id_counter = -1 self.max_steps = max_steps self.grid.add_nest(nest_position) - self.num_max_agents = num_max_agents - self.num_new_recruits = 5 - self.decay_rates : dict[str, float] = {"A" :0.01, - "B": 0.01, - } - - for agent_id in self.get_unique_ids(num_initial_roamers): - if self.schedule.get_agent_count() < self.num_max_agents: + for agent_id in self.get_unique_ids(N_0): + if self.schedule.get_agent_count() < self.N_m: agent = RandomWalkerAnt(unique_id=agent_id, model=self, look_for_pheromone="A", drop_pheromone="A") self.schedule.add(agent) self.grid.place_agent(agent, pos=nest_position) - for _ in range(num_food_sources): + for _ in range(N_f): self.grid.add_food(food_size) self.datacollector = DataCollector( - model_reporters={}, + # model_reporters={"agent_dens": lambda m: m.agent_density()}, + model_reporters = {"pheromone_a": lambda m: m.grid.fields["A"], + "pheromone_b": lambda m: m.grid.fields["B"], + }, agent_reporters={} ) self.datacollector.collect(self) # keep at end of __init___ @@ -68,7 +159,7 @@ class ActiveWalkerModel(Model): # apply decay rate on pheromone levels for key in ("A", "B"): field = self.grid.fields[key] - self.grid.fields[key] = field - self.decay_rates[key]*field + self.grid.fields[key] = field - self.gamma*field self.datacollector.collect(self) diff --git a/multihex.py b/multihex.py index e1ebae2..5f80d5b 100644 --- a/multihex.py +++ b/multihex.py @@ -107,7 +107,8 @@ class MultiHexGridScalarFields(MultiHexGrid): def is_food(self, pos): assert('food' in self.fields.keys()) - return bool(self.fields['food'][pos]) + # account for potential float imprecision and use epsilon = 1e-3 + return self.fields['food'][pos] > 1e-3 def add_food(self, size : int , pos=None): """ @@ -127,7 +128,7 @@ class MultiHexGridScalarFields(MultiHexGrid): while(self.is_nest(pos) or self.is_food(pos)): pos = select_random_place() - self.fields['food'][pos] = size + self.fields['food'][pos] = int(size) def is_nest(self, pos : Coordinate) -> bool: assert('nests' in self.fields.keys()) diff --git a/server.py b/server.py index 2a2a159..ad5924c 100755 --- a/server.py +++ b/server.py @@ -22,10 +22,12 @@ def setup(params=None): # Set the model parameters if params is None: params = { + "max_steps": 3000, "width": 50, "height": 50, - "num_max_agents" : 100, + "N_m" : 100, "nest_position" : (25,25), - "num_initial_roamers" : 5, + "N_0" : 5, + "resistance_map_type": "perlin", } @@ -85,12 +87,9 @@ def setup(params=None): "Color": col, } - def get_max_grid_val(model, key): - return np.max(model.grid.fields[key]) - - def portray_pheromone_density(model, pos, norm): - col_a = get_color(level=model.grid.fields["A"][pos], normalization=norm) - col_b = get_color(level=model.grid.fields["B"][pos], normalization=norm) + def portray_resistance_map(model, pos, norm=1): + col = get_color(level=model.grid.fields['res'][pos], normalization=norm) + col = f"rgb({col}, {col}, {col})" return { "Shape": "hex", "r": 1, @@ -98,7 +97,26 @@ def setup(params=None): "Layer": 0, "x": pos[0], "y": pos[1], - "Color": f"rgb({col_a}, {col_b}, 255)" + "Color": col, + } + + def get_max_grid_val(model, key): + return np.max(model.grid.fields[key]) + + def portray_pheromone_density(model, pos, norm): + col_a = get_color(level=model.grid.fields["A"][pos], normalization=norm) + col_b = get_color(level=model.grid.fields["B"][pos], normalization=norm) + res_min, res_max = np.min(model.grid.fields['res']), np.max(model.grid.fields['res']) + ease = 1 - model.grid.fields['res'][pos] + col_ease = get_color(level=ease, normalization=np.max(model.grid.fields['res'])) + return { + "Shape": "hex", + "r": 1, + "Filled": "true", + "Layer": 0, + "x": pos[0], + "y": pos[1], + "Color": f"rgb({col_a}, {col_b}, {col_ease})" } @@ -109,6 +127,9 @@ def setup(params=None): grid_ants = CanvasHexGridMultiAgents(portray_ant_density, width, height, width*pixel_ratio, height*pixel_ratio, norm_method=lambda m: 5) + grid_resistance_map = CanvasHexGridMultiAgents(portray_resistance_map, + width, height, width*pixel_ratio, height*pixel_ratio, + norm_method=lambda m: 1) def norm_ants(model): return 5 @@ -127,14 +148,19 @@ def setup(params=None): [lambda m: "

Ant density

Nest: Red, Food: Green
", grid_ants, lambda m: f"
Normalization Value: {norm_ants(m)}
", - lambda m: "

Pheromone Density

Pheromone A: Cyan, Pheromone B: Magenta
", + lambda m: "

Pheromone Density

Pheromone A: Cyan, Pheromone B: Magenta, Resistance Map: Yellow
", grid_pheromones, - lambda m: f"
Normalization Value: {norm_pheromones(m)}
" + lambda m: f"
Normalization Value: {norm_pheromones(m)}
", ], "Active Random Walker Ants", params) if __name__ == "__main__": - server = setup() + from model import kwargs_paper_setup1 + kwargs_paper1_perlin = kwargs_paper_setup1 + kwargs_paper1_perlin["height"] = 50 + kwargs_paper1_perlin["width"] = 50 + kwargs_paper1_perlin["resistance_map_type"] = "perlin" + server = setup(params=kwargs_paper1_perlin) server.launch() """