Implement El Farol model

Co-authored-by: Daniel Xu <[email protected]>

examples/el_farol/README.md

@@ -0,0 +1,22 @@
+# El Farol
+
+This folder contains an implementation of El Farol restaurant model. Agents (restaurant customers) decide whether to go to the restaurant or not based on their memory and reward from previous trials. Implications from the model have been used to explain how individual decision-making affects overall performance and fluctuation.
+
+The implementation is based on Fogel 1999 (in particular the calculation of the prediction), which is a refinement over Arthur 1994.
+
+## How to Run
+
+Launch the model: You can run the model and perform analysis in el_farol.ipynb.
+You can test the model itself by running `pytest tests.py`.
+
+## Files
+* [el_farol.ipynb](el_farol.ipynb): Run the model and visualization in a Jupyter notebook
+* [el_farol/model.py](el_farol/model.py): Core model file.
+* [el_farol/agents.py](el_farol/agents.py): The agent class.
+* [tests.py](tests.py): Tests to ensure the model is consistent with Arthur 1994, Fogel 1996.
+
+## Further Reading
+
+1. W. Brian Arthur Inductive Reasoning and Bounded Rationality (1994) https://www.jstor.org/stable/2117868
+1. D.B. Fogel, K. Chellapilla, P.J. Angeline Inductive reasoning and bounded rationality reconsidered (1999)
+1. NetLogo implementation of the El Farol bar problem https://ccl.northwestern.edu/netlogo/models/ElFarol

examples/el_farol/el_farol.ipynb

@@ -0,0 +1,157 @@
+{
+ "cells": [
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import matplotlib.pyplot as plt\n",
+    "import numpy as np\n",
+    "import seaborn as sns\n",
+    "\n",
+    "from el_farol.model import ElFarolBar"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "memory_sizes = [5, 10, 20]\n",
+    "crowd_threshold = 60\n",
+    "models = [\n",
+    "    ElFarolBar(N=100, crowd_threshold=crowd_threshold, memory_size=m)\n",
+    "    for m in memory_sizes\n",
+    "]\n",
+    "for model in models:\n",
+    "    for i in range(100):\n",
+    "        model.step()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# You should observe that the attendance converges to 60.\n",
+    "_, axs = plt.subplots(1, 3, figsize=(10, 3))\n",
+    "for idx, model in enumerate(models):\n",
+    "    ax = axs[idx]\n",
+    "    plt.sca(ax)\n",
+    "    df = model.datacollector.get_model_vars_dataframe()\n",
+    "    sns.lineplot(data=df, x=df.index, y=\"Customers\", ax=ax)\n",
+    "    ax.set(\n",
+    "        xlabel=\"Step\",\n",
+    "        ylabel=\"Attendance\",\n",
+    "        title=f\"Memory size = {memory_sizes[idx]}\",\n",
+    "        ylim=(20, 80),\n",
+    "    )\n",
+    "    plt.axhline(crowd_threshold, color=\"tab:red\")\n",
+    "    plt.tight_layout()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "for idx, memory_size in enumerate(memory_sizes):\n",
+    "    model = models[idx]\n",
+    "    df = model.datacollector.get_agent_vars_dataframe()\n",
+    "    sns.lineplot(\n",
+    "        x=df.index.levels[0],\n",
+    "        y=df.Utility.groupby(\"Step\").mean(),\n",
+    "        label=str(memory_size),\n",
+    "    )\n",
+    "plt.legend(title=\"Memory size\");"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Decisions made on across trials\n",
+    "fix, axs = plt.subplots(1, 3, figsize=(12, 4))\n",
+    "for idx, memory_size in enumerate(memory_sizes):\n",
+    "    plt.sca(axs[idx])\n",
+    "    df = models[idx].datacollector.get_agent_vars_dataframe()\n",
+    "    df.reset_index(inplace=True)\n",
+    "    ax = sns.heatmap(df.pivot(index=\"AgentID\", columns=\"Step\", values=\"Attendance\"))\n",
+    "    ax.set(title=f\"Memory size = {memory_size}\")\n",
+    "    plt.tight_layout()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Next, we experiment with varying the number of strategies\n",
+    "num_strategies_list = [5, 10, 20]\n",
+    "crowd_threshold = 60\n",
+    "models = [\n",
+    "    ElFarolBar(N=100, crowd_threshold=crowd_threshold, num_strategies=ns)\n",
+    "    for ns in num_strategies_list\n",
+    "]\n",
+    "for model in models:\n",
+    "    for i in range(100):\n",
+    "        model.step()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Attendance of the bar based on the number of strategies\n",
+    "_, axs = plt.subplots(1, 3, figsize=(10, 3))\n",
+    "for idx, num_strategies in enumerate(num_strategies_list):\n",
+    "    model = models[idx]\n",
+    "    ax = axs[idx]\n",
+    "    plt.sca(ax)\n",
+    "    df = model.datacollector.get_model_vars_dataframe()\n",
+    "    sns.lineplot(data=df, x=df.index, y=\"Customers\", ax=ax)\n",
+    "    ax.set(\n",
+    "        xlabel=\"Trial\",\n",
+    "        ylabel=\"Attendance\",\n",
+    "        title=f\"Number of Strategies = {num_strategies}\",\n",
+    "        ylim=(20, 80),\n",
+    "    )\n",
+    "    plt.axhline(crowd_threshold, color=\"tab:red\")\n",
+    "    plt.tight_layout()"
+   ]
+  }
+ ],
+ "metadata": {
+  "interpreter": {
+   "hash": "18b8a6ab22c23ac88fce14986952a46f0d293914064547c699eac09fb58cfe0f"
+  },
+  "kernelspec": {
+   "display_name": "Python 3 (ipykernel)",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.11.6"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 4
+}

examples/el_farol/el_farol/agents.py

@@ -0,0 +1,63 @@
+import mesa
+import numpy as np
+
+
+class BarCustomer(mesa.Agent):
+    def __init__(self, unique_id, model, memory_size, crowd_threshold, num_strategies):
+        super().__init__(unique_id, model)
+        # Random values from -1.0 to 1.0
+        self.strategies = np.random.rand(num_strategies, memory_size + 1) * 2 - 1
+        self.best_strategy = self.strategies[0]
+        self.attend = False
+        self.memory_size = memory_size
+        self.crowd_threshold = crowd_threshold
+        self.utility = 0
+        self.update_strategies()
+
+    def step(self):
+        prediction = self.predict_attendance(
+            self.best_strategy, self.model.history[-self.memory_size :]
+        )
+        if prediction <= self.crowd_threshold:
+            self.attend = True
+            self.model.attendance += 1
+        else:
+            self.attend = False
+
+    def update_strategies(self):
+        # Pick the best strategy based on new history window
+        best_score = float("inf")
+        for strategy in self.strategies:
+            score = 0
+            for week in range(self.memory_size):
+                last = week + self.memory_size
+                prediction = self.predict_attendance(
+                    strategy, self.model.history[week:last]
+                )
+                score += abs(self.model.history[last] - prediction)
+            if score <= best_score:
+                best_score = score
+                self.best_strategy = strategy
+        should_attend = self.model.history[-1] <= self.crowd_threshold
+        if should_attend != self.attend:
+            self.utility -= 1
+        else:
+            self.utility += 1
+
+    def predict_attendance(self, strategy, subhistory):
+        # This is extracted from the source code of the model in
+        # https://ccl.northwestern.edu/netlogo/models/ElFarol.
+        # This reports an agent's prediction of the current attendance
+        # using a particular strategy and portion of the attendance history.
+        # More specifically, the strategy is then described by the formula
+        # p(t) = x(t - 1) * a(t - 1) + x(t - 2) * a(t - 2) +..
+        #      ... + x(t - memory_size) * a(t - memory_size) + c * 100,
+        # where p(t) is the prediction at time t, x(t) is the attendance of the
+        # bar at time t, a(t) is the weight for time t, c is a constant, and
+        # MEMORY-SIZE is an external parameter.
+
+        # The first element of the strategy is the constant, c, in the
+        # prediction formula. one can think of it as the the agent's prediction
+        # of the bar's attendance in the absence of any other data then we
+        # multiply each week in the history by its respective weight.
+        return strategy[0] * 100 + np.dot(strategy[1:], subhistory)

examples/el_farol/el_farol/model.py

@@ -0,0 +1,45 @@
+import mesa
+import numpy as np
+
+from .agents import BarCustomer
+
+
+class ElFarolBar(mesa.Model):
+    def __init__(
+        self,
+        crowd_threshold=60,
+        num_strategies=10,
+        memory_size=10,
+        width=100,
+        height=100,
+        N=100,
+    ):
+        self.running = True
+        self.num_agents = N
+        self.schedule = mesa.time.RandomActivation(self)
+
+        # Initialize the previous attendance randomly so the agents have a history
+        # to work with from the start.
+        # The history is twice the memory, because we need at least a memory
+        # worth of history for each point in memory to test how well the
+        # strategies would have worked.
+        self.history = np.random.randint(0, 100, size=memory_size * 2).tolist()
+        self.attendance = self.history[-1]
+        for i in range(self.num_agents):
+            a = BarCustomer(i, self, memory_size, crowd_threshold, num_strategies)
+            self.schedule.add(a)
+        self.datacollector = mesa.DataCollector(
+            model_reporters={"Customers": "attendance"},
+            agent_reporters={"Utility": "utility", "Attendance": "attend"},
+        )
+
+    def step(self):
+        self.datacollector.collect(self)
+        self.attendance = 0
+        self.schedule.step()
+        # We ensure that the length of history is constant
+        # after each step.
+        self.history.pop(0)
+        self.history.append(self.attendance)
+        for agent in self.schedule.agents:
+            agent.update_strategies()

examples/el_farol/requirements.txt

@@ -0,0 +1,5 @@
+jupyter
+matplotlib
+mesa
+numpy
+seaborn

examples/el_farol/tests.py

@@ -0,0 +1,19 @@
+import numpy as np
+from el_farol.model import ElFarolBar
+
+np.random.seed(1)
+crowd_threshold = 60
+
+
+def test_convergence():
+    # Testing that the attendance converges to crowd_threshold
+    attendances = []
+    for _ in range(10):
+        model = ElFarolBar(N=100, crowd_threshold=crowd_threshold, memory_size=10)
+        for _ in range(100):
+            model.step()
+        attendances.append(model.attendance)
+    mean = np.mean(attendances)
+    standard_deviation = np.std(attendances)
+    deviation = abs(mean - crowd_threshold)
+    assert deviation < standard_deviation