Implement El Farol model

Co-authored-by: Daniel Xu <tx29@uw.edu>

examples/el_farol/README.md

@@ -0,0 +1,29 @@
+# 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.
+
+This example has 2 versions of the model: the pure Mesa version without IBL (named as ElFarolBar), and the version with IBL that integrates with PyIBL (named as ElFarolBarIBLT). The latter version replicates the result of Kumar 2016. The base version without IBL is based on Fogel 1999 (in particular the calculation of the prediction), which is a refinement over Arthur 1994.
+
+The IBL version of the model demonstrates how to deploy a cognitive model (Instance-Based Learning) under the Mesa environment. IBL model reflects the recency and frequency effect in decision-making with memory. Agent actively learns from the environment and updates their preference(blending value) for each decision. IBL model could be used as a substitute for an agent whose decision-making is more realistic and closer to human decision-making.
+TODO: The first plot in el_farol_iblt.ipynb does not match figure 1 in Kumar 2016.
+
+
+## How to Run
+
+Launch the model: Please check el_farol.ipynb for more information.
+Please see this [link](http://pyibl.ddmlab.com/) to install pyibl package.
+
+## Files
+* [el_farol.ipynb](el_farol.ipynb): Run the model and visualization in a Jupyter notebook
+* [el_farol_iblt.ipynb](el_farol_iblt.ipynb): Run the IBLT 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 and also contain a cognitive model for el_farol problem.
+* [tests.py](tests.py): Tests to ensure the model is consistent with Arthur 1994, Fogel 1996, Kumar 2016.
+
+## Further Reading
+
+=======
+[1] W. Brian Arthur Inductive Reasoning and Bounded Rationality (1994) https://www.jstor.org/stable/2117868
+[2] D.B. Fogel, K. Chellapilla, P.J. Angeline Inductive reasoning and bounded rationality reconsidered (1999)
+[3] NetLogo implementation of the El Farol bar problem https://ccl.northwestern.edu/netlogo/models/ElFarol
+[3] Kumar, Shikhar, and Cleotilde Gonzalez. "Heterogeneity of Memory Decay and Collective Learning in the El Farol Bar Problem." (2016).

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,109 @@
+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 + sum(strategy[1:] * subhistory)
+
+
+class BarCustomerIBLT(mesa.Agent):
+    """
+    This is BarCustomer but implemented using PyIBL
+    """
+
+    def __init__(self, unique_id, model, decay, crowd_threshold):
+        super().__init__(unique_id, model)
+
+        import pyibl
+
+        self.agent = pyibl.Agent(
+            name="BarCustomer",
+            attributes=["Attendance"],
+            decay=decay,
+            noise=np.random.uniform(0.2, 0.8),
+        )
+        self.agent.default_utility = 10
+        self.utility = 0
+        self.decay = decay
+        self.crowd_threshold = crowd_threshold
+        # The step() at initialization is necessary because the agent respond
+        # needs the choose method to be executed beforehand.
+        self.step()
+        self.update_strategies()
+
+    def step(self):
+        choice = self.agent.choose(["Attend", "Not Attend"])
+        if choice == "Attend":
+            self.attend = True
+            self.model.attendance += 1
+        else:
+            self.attend = False
+
+    def update_strategies(self):
+        """
+        Update blending value for IBL agent. the if statement is the same as the if statement in BarCustomer agent
+        """
+        should_attend = self.model.history[-1] <= self.crowd_threshold
+        if should_attend != self.attend:
+            self.agent.respond(-1)
+            self.utility -= 1
+        else:
+            self.agent.respond(1)
+            self.utility += 1

examples/el_farol/el_farol/model.py

@@ -0,0 +1,78 @@
+import mesa
+import numpy as np
+
+from .agents import BarCustomer, BarCustomerIBLT
+
+
+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()
+
+
+class ElFarolBarIBLT(ElFarolBar):
+    def __init__(
+        self,
+        crowd_threshold=60,
+        decay_portions=None,
+        memory_size=10,
+        width=100,
+        height=100,
+        N=100,
+    ):
+        self.running = True
+        self.num_agents = N
+        self.schedule = mesa.time.RandomActivation(self)
+        self.history = np.random.randint(0, 100, size=memory_size * 2).tolist()
+        self.attendance = self.history[0]
+        if decay_portions is None:
+            decay_portions = {1: 1}
+        i = 0
+        for decay, portion in decay_portions.items():
+            for _ in range(int(self.num_agents * portion)):
+                a = BarCustomerIBLT(i, self, decay, crowd_threshold)
+                self.schedule.add(a)
+                i += 1
+        self.datacollector = mesa.DataCollector(
+            model_reporters={"Customers": "attendance"},
+            agent_reporters={
+                "Utility": "utility",
+                "Decay": "decay",
+                "Attendance": "attend",
+            },
+        )