General Debugging. Moved to nomenclature from paper for ease of

understanding
Resistance map implementation

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()

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)

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()
 
 
 

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)
 

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())

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: "<h3>Ant density</h3><h5>Nest: Red, Food: Green</h5>",
                           grid_ants,
                           lambda m: f"<h5>Normalization Value: {norm_ants(m)}</h5>",
-                          lambda m: "<h3>Pheromone Density</h3><h5>Pheromone A: Cyan, Pheromone B: Magenta</h5>",
+                          lambda m: "<h3>Pheromone Density</h3><h5>Pheromone A: Cyan, Pheromone B: Magenta, Resistance Map: Yellow</h5>",
                           grid_pheromones,
-                          lambda m: f"<h5>Normalization Value: {norm_pheromones(m)}</h5>"
+                          lambda m: f"<h5>Normalization Value: {norm_pheromones(m)}</h5>",
                           ],
                            "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()
 
 """