{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Dynamical Systems Models: The Logistic ODE and Behavioral Acquisition\n",
    "\n",
    "In this lab you will work with the logistic ordinary differential equation (ODE) as a model of behavioral acquisition:\n",
    "\n",
    "$$\\frac{dx}{dt} = r \\cdot x \\cdot \\left(1 - \\frac{x}{K}\\right)$$\n",
    "\n",
    "where:\n",
    "- $x(t)$ is the response rate at time $t$\n",
    "- $r$ is the intrinsic growth rate (how quickly behavior accelerates early in acquisition)\n",
    "- $K$ is the carrying capacity (the asymptotic response rate)\n",
    "\n",
    "**Learning objectives:**\n",
    "- Implement numerical integration using Euler's method\n",
    "- Analyze equilibrium points and their stability\n",
    "- Construct phase portraits\n",
    "- Fit dynamical models to empirical data\n",
    "\n",
    "**Required packages:** `numpy`, `matplotlib`, `scipy`"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Setup\n",
    "\n",
    "Run the cell below to import the libraries you will need."
   ]
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy.optimize import curve_fit\n",
    "\n",
    "plt.rcParams['figure.figsize'] = (10, 6)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Task 1: Define the Logistic ODE and Implement Euler's Method\n",
    "\n",
    "First, define a Python function for the right-hand side of the logistic ODE:\n",
    "\n",
    "$$f(x) = r \\cdot x \\cdot \\left(1 - \\frac{x}{K}\\right)$$\n",
    "\n",
    "Then implement Euler's method to solve the ODE numerically. Recall that Euler's method updates the state as:\n",
    "\n",
    "$$x_{n+1} = x_n + \\Delta t \\cdot f(x_n)$$\n",
    "\n",
    "Write a function `euler_logistic(r, K, x0, dt, n_steps)` that returns arrays of time values and corresponding $x$ values.\n",
    "\n",
    "Use the following parameters for your first run: $r = 0.2$, $K = 45$, $x_0 = 1.5$, $\\Delta t = 0.1$, and enough steps to cover 30 time units."
   ]
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "# Define logistic_ode(x, r, K) and euler_logistic(r, K, x0, dt, n_steps)\n"
   ]
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "# Run Euler's method (r=0.2, K=45, x0=1.5, dt=0.1, 30 time units) and plot\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Task 2: Explore the Effect of Parameters and Initial Conditions\n",
    "\n",
    "Create a figure with two subplots:\n",
    "\n",
    "**Subplot 1 - Varying initial conditions:** Fix $r = 0.2$ and $K = 45$. Plot the acquisition curve for $x_0 \\in \\{0.5, 1.5, 5.0, 15.0\\}$. All trajectories should converge to $K$, but they will take different amounts of time to get there.\n",
    "\n",
    "**Subplot 2 - Varying growth rate:** Fix $K = 45$ and $x_0 = 1.5$. Plot the acquisition curve for $r \\in \\{0.1, 0.2, 0.4, 0.8\\}$. How does $r$ affect the steepness and timing of the inflection point?\n",
    "\n",
    "Label your axes and include legends."
   ]
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "# Two subplots: vary x0, then vary r\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Question:** In behavioral terms, what does $r$ represent? What does $K$ represent? How might these parameters relate to schedule of reinforcement and response ceiling effects?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*Your answer here.*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Task 3: Find Equilibrium Points Analytically\n",
    "\n",
    "The equilibrium points (also called fixed points or steady states) are values of $x^*$ where $\\frac{dx}{dt} = 0$.\n",
    "\n",
    "Set $f(x^*) = r \\cdot x^* \\cdot (1 - x^*/K) = 0$ and solve for $x^*$.\n",
    "\n",
    "1. Write out the equilibrium points analytically in the markdown cell below.\n",
    "2. Verify your answer by evaluating `logistic_ode(x_star, r, K)` for each equilibrium and confirming the result is zero (or numerically very close to zero)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Equilibrium points:**\n",
    "\n",
    "*Your analytical derivation here.*"
   ]
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "# Verify f(x*) = 0 at each equilibrium\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Task 4: Linear Stability Analysis\n",
    "\n",
    "To determine whether each equilibrium is stable or unstable, compute the derivative of $f(x)$ with respect to $x$:\n",
    "\n",
    "$$f'(x) = \\frac{d}{dx}\\left[r \\cdot x \\cdot \\left(1 - \\frac{x}{K}\\right)\\right]$$\n",
    "\n",
    "Evaluate $f'(x^*)$ at each equilibrium point:\n",
    "- If $f'(x^*) < 0$, the equilibrium is **stable** (perturbations decay).\n",
    "- If $f'(x^*) > 0$, the equilibrium is **unstable** (perturbations grow).\n",
    "\n",
    "1. Derive $f'(x)$ analytically.\n",
    "2. Evaluate it at each equilibrium for $r = 0.2$ and $K = 45$.\n",
    "3. State which equilibrium is stable and which is unstable, and explain what this means for the acquisition of behavior."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Derivation of $f'(x)$:**\n",
    "\n",
    "*Your work here.*"
   ]
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "# Define f'(x) and evaluate at each equilibrium\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Task 5: Create a Phase Portrait\n",
    "\n",
    "A phase portrait plots $\\frac{dx}{dt}$ on the y-axis against $x$ on the x-axis. This shows at a glance how the rate of change depends on the current state.\n",
    "\n",
    "1. Create an array of $x$ values from $-2$ to $55$ (going slightly beyond $K$).\n",
    "2. Compute $f(x)$ for each value.\n",
    "3. Plot $f(x)$ vs. $x$.\n",
    "4. Mark the equilibrium points with dots.\n",
    "5. Add arrows along the x-axis indicating the direction of flow:\n",
    "   - Where $f(x) > 0$, $x$ is increasing (rightward arrows).\n",
    "   - Where $f(x) < 0$, $x$ is decreasing (leftward arrows).\n",
    "6. Draw a horizontal dashed line at $f(x) = 0$.\n",
    "\n",
    "**Hint:** Use `plt.annotate` with `arrowprops` or `plt.quiver` for the arrows."
   ]
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "# Phase portrait: f(x) vs x, mark equilibria, add flow arrows\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Question:** Looking at the phase portrait, explain in your own words why $x^* = K$ is a stable equilibrium and $x^* = 0$ is unstable (for $x > 0$). What happens if $x$ is slightly above $K$?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*Your answer here.*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Task 6: The Analytical Solution\n",
    "\n",
    "The logistic ODE has a well-known analytical (closed-form) solution:\n",
    "\n",
    "$$x(t) = \\frac{K}{1 + \\left(\\frac{K - x_0}{x_0}\\right) e^{-rt}}$$\n",
    "\n",
    "1. Implement this as a Python function.\n",
    "2. Plot it alongside the Euler numerical solution from Task 1 to verify they agree.\n",
    "3. Try a larger time step (e.g., $\\Delta t = 1.0$) in Euler's method and compare. How does the step size affect accuracy?"
   ]
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "# Define logistic_analytical(t, r, K, x0)\n"
   ]
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "# Compare analytical solution with Euler at dt=0.1 and dt=1.0\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Task 7: Fit the Logistic Model to Empirical Acquisition Data\n",
    "\n",
    "Load `acquisition_data.csv`, which contains response rates across 30 sessions for a single participant learning a new operant response.\n",
    "\n",
    "1. Load the data and plot the raw acquisition curve.\n",
    "2. Use `scipy.optimize.curve_fit` to fit the analytical logistic solution to the data. You will need to define a wrapper function that takes `t` as the first argument and the parameters to fit ($r$, $K$, $x_0$) as subsequent arguments.\n",
    "3. Provide reasonable initial guesses: $r_0 = 0.3$, $K_0 = 40$, $x_{0,0} = 1.0$.\n",
    "4. Report the fitted parameter values and their standard errors.\n",
    "5. Plot the data and the fitted curve together."
   ]
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "acq_data = pd.read_csv('acquisition_data.csv')\n",
    "# Plot the raw acquisition curve\n"
   ]
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "# Define the curve_fit wrapper and fit (p0 = [0.3, 40, 1.0])\n"
   ]
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "# Report fitted parameters (+/- std errors) and plot data vs fit\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Task 8: Assess Model Fit\n",
    "\n",
    "1. Compute the residuals (observed - predicted) and plot them against session number.\n",
    "2. Compute $R^2$ for the fit.\n",
    "3. Are the residuals randomly scattered, or do you see systematic patterns? What might this suggest about the adequacy of the logistic model for these data?"
   ]
  },
  {
   "cell_type": "code",
   "metadata": {},
   "execution_count": null,
   "outputs": [],
   "source": [
    "# Residuals vs session, and R-squared\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Task 9: Behavioral Interpretation\n",
    "\n",
    "Discuss the following in the cell below:\n",
    "\n",
    "1. What is the estimated carrying capacity $K$? How does this compare to what you would expect as a maximum response rate for an operant response?\n",
    "2. What is the estimated growth rate $r$? How might this parameter change under different reinforcement schedules (e.g., continuous vs. intermittent reinforcement)?\n",
    "3. What is the estimated initial response rate $x_0$? Does it make sense as a starting level of behavior?\n",
    "4. What are the limitations of using the logistic model for behavioral acquisition? Can you think of behavioral phenomena it would fail to capture (e.g., extinction, ratio strain, contrast effects)?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*Your discussion here.*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Reflection\n",
    "\n",
    "In a few sentences, address the following:\n",
    "\n",
    "1. How does thinking about behavior in terms of differential equations change your perspective compared to static curve-fitting?\n",
    "2. What is the advantage of having both a numerical method (Euler) and an analytical solution? When might you only have one or the other?\n",
    "3. How could you extend the logistic model to capture more complex behavioral dynamics (e.g., extinction following acquisition, or competing responses)?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*Your reflection here.*"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "name": "python"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}