{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Pynamical: demo of the logistic map and bifurcation diagrams\n",
    "\n",
    "**Citation info**: Boeing, G. 2016. \"[Visual Analysis of Nonlinear Dynamical Systems: Chaos, Fractals, Self-Similarity and the Limits of Prediction](http://geoffboeing.com/publications/nonlinear-chaos-fractals-prediction/).\" *Systems*, 4 (4), 37. doi:10.3390/systems4040037.\n",
    "\n",
    "Pynamical documentation: http://pynamical.readthedocs.org\n",
    "\n",
    "This notebook implements a logistic map and plots its results, bifurcation diagrams, and phase diagrams"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import IPython.display as IPdisplay\n",
    "import matplotlib.cm as cm\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import pynamical\n",
    "from pynamical import simulate, bifurcation_plot, save_fig\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "title_font = pynamical.get_title_font()\n",
    "label_font = pynamical.get_label_font()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## First, let's see the population values the logistic map produces for a range of growth rate parameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# run the logistic model for 20 generations for 7 growth rates between 0.5 and 3.5 then view the output\n",
    "pops = simulate(num_gens=20, rate_min=0.5, rate_max=3.5, num_rates=7)\n",
    "pops.applymap(lambda x: '{:03.3f}'.format(x))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_colors(cmap, n, start=0., stop=1., alpha=1., reverse=False):\n",
    "    '''return n-length list of rgba colors from the passed colormap name and alpha,\n",
    "       limit extent by start/stop values and reverse list order if flag is true'''\n",
    "    colors = [cm.get_cmap(cmap)(x) for x in np.linspace(start, stop, n)]\n",
    "    colors = [(r, g, b, alpha) for r, g, b, _ in colors]\n",
    "    return list(reversed(colors)) if reverse else colors"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# plot the results of the logistic map run for these 7 different growth rates\n",
    "#color_list = ['#cc00cc', '#4B0082', '#0066cc', '#33cc00', '#cccc33', '#ff9900', '#ff0000']\n",
    "color_list = get_colors('viridis', n=len(pops.columns), start=0., stop=1)\n",
    "for color, rate in reversed(list(zip(color_list, pops.columns))):\n",
    "    ax = pops[rate].plot(kind='line', figsize=[10, 6], linewidth=2.5, alpha=0.95, c=color)\n",
    "ax.grid(True)\n",
    "ax.set_ylim([0, 1])\n",
    "ax.legend(title='Growth Rate', loc=3, bbox_to_anchor=(1, 0.525))\n",
    "ax.set_title('Logistic Model Results by Growth Rate', fontproperties=title_font)\n",
    "ax.set_xlabel('Generation', fontproperties=label_font)\n",
    "ax.set_ylabel('Population', fontproperties=label_font)\n",
    "\n",
    "save_fig('logistic-map-growth-rates')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Now let's visualize the system attractors for a large range of growth rate parameters, using bifurcation diagrams"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# run the model for 100 generations across 1000 growth rate steps from 0 to 4 then plot the bifurcation diagram\n",
    "pops = simulate(num_gens=100, rate_min=0, rate_max=4, num_rates=1000, num_discard=1)\n",
    "bifurcation_plot(pops, filename='logistic-map-bifurcation-0')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# plot the bifurcation diagram for 200 generations, but this time throw out the first 100 rows\n",
    "# 200-100=100, so we still have 100 generations in the plot, just like in the previous cell\n",
    "# this will show us only the attractors (aka, the values that each growth rate settles on over time)\n",
    "pops = simulate(num_gens=100, rate_min=0, rate_max=4, num_rates=1000, num_discard=100)\n",
    "bifurcation_plot(pops, filename='logistic-map-bifurcation-1')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# run the model for 300 generations across 1,000 growth rate steps from 2.8 to 4, and plot the bifurcation diagram\n",
    "# this plot is a zoomed-in look at the first plot and shows the period-doubling path to chaos\n",
    "pops = simulate(num_gens=100, rate_min=2.8, rate_max=4, num_rates=1000, num_discard=200, initial_pop=0.1)\n",
    "bifurcation_plot(pops, xmin=2.8, xmax=4, filename='logistic-map-bifurcation-2')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# run the model for 200 generations across 1,000 growth rate steps from 3.7 to 3.9, and plot the bifurcation diagram\n",
    "# this plot is a zoomed-in look at the first plot and shows more detail in the chaotic regimes\n",
    "pops = simulate(num_gens=100, rate_min=3.7, rate_max=3.9, num_rates=1000, num_discard=100)\n",
    "bifurcation_plot(pops, xmin=3.7, xmax=3.9, filename='logistic-map-bifurcation-3')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## In the chaotic regime (r=3.6 to 4=4.0), the system has a strange attractor with fractal structure"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# run the model for 500 generations across 1,000 growth rate steps from 3.84 to 3.856, and plot the bifurcation diagram\n",
    "# throw out the first 300 generations, so we end up with 200 generations in the plot\n",
    "# this plot is a zoomed-in look at the first plot and shows the same structure we saw at the macro-level\n",
    "pops = simulate(num_gens=200, rate_min=3.84, rate_max=3.856, num_rates=1000, num_discard=300)\n",
    "bifurcation_plot(pops, xmin=3.84, xmax=3.856, ymin=0.445, ymax=0.552, filename='logistic-map-bifurcation-4')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Now let's visualize the system's sensitive dependence on initial conditions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# plot the numeric output of the logistic model for growth rates of 3.9 and 3.90001\n",
    "# this demonstrates sensitive dependence on the parameter\n",
    "rate1 = 3.9\n",
    "rate2 = rate1 + 0.00001\n",
    "\n",
    "pops = simulate(num_gens=40, rate_min=rate1, rate_max=rate2, num_rates=2)\n",
    "ax = pops.plot(kind='line', figsize=[10, 6], linewidth=3, alpha=0.6, style=['#003399','#cc0000'])\n",
    "ax.grid(True)\n",
    "ax.set_title('Logistic Model Results by Growth Rate', fontproperties=title_font)\n",
    "ax.set_xlabel('Generation', fontproperties=label_font)\n",
    "ax.set_ylabel('Population', fontproperties=label_font)\n",
    "ax.legend(title='Growth Rate', loc=3)\n",
    "\n",
    "save_fig('logistic-map-parameter-sensitivity')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# plot the numeric output of the logistic model at growth rate 3.9 for 2 similar starting population values\n",
    "# this demonstrates sensitive dependence on initial conditions, as they diverge through chaos\n",
    "r = 3.9\n",
    "pops1 = simulate(num_gens=55, rate_min=r, rate_max=4.0, num_rates=1, initial_pop=0.5)\n",
    "pops2 = simulate(num_gens=55, rate_min=r, rate_max=4.0, num_rates=1, initial_pop=0.50001)\n",
    "pops = pd.concat([pops1, pops2], axis=1)\n",
    "pops.columns = ['0.5', '0.50001']\n",
    "ax = pops.plot(kind='line', figsize=[10, 6], linewidth=3, alpha=0.6, style=['#003399','#cc0000'])\n",
    "ax.grid(True)\n",
    "ax.set_title('Logistic Model Results by Initial Conditions, r={}'.format(r), fontproperties=title_font)\n",
    "ax.set_xlabel('Generation', fontproperties=label_font)\n",
    "ax.set_ylabel('Population', fontproperties=label_font)\n",
    "ax.legend(title='Initial Population', loc=3)\n",
    "\n",
    "save_fig('logistic-map-initial-conditions')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# plot the numeric output of the logistic model at growth rate 3.65 for 2 similar starting population values\n",
    "# this demonstrates how very similar conditions do not diverge when the rate is not chaotic\n",
    "r = 3.65\n",
    "pops1 = simulate(num_gens=55, rate_min=r, num_rates=1, initial_pop=0.5)\n",
    "pops2 = simulate(num_gens=55, rate_min=r, num_rates=1, initial_pop=0.50001)\n",
    "pops = pd.concat([pops1, pops2], axis=1)\n",
    "pops.columns = ['0.5', '0.50001']\n",
    "ax = pops.plot(kind='line', figsize=[10, 6], linewidth=3, alpha=0.6, style=['#003399','#cc0000'])\n",
    "ax.grid(True)\n",
    "ax.set_title('Logistic Model Results by Initial Conditions, r={}'.format(r), fontproperties=title_font)\n",
    "ax.set_xlabel('Generation', fontproperties=label_font)\n",
    "ax.set_ylabel('Population', fontproperties=label_font)\n",
    "ax.legend(title='Initial Population', loc=3)\n",
    "\n",
    "save_fig('logistic-map-initial-conditions-stable')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## In part 2, I look at phase diagrams that let us visualize our strange attractors and disambiguate chaos from random noise:\n",
    "\n",
    "[pynamical-demo-phase-diagrams.ipynb](pynamical-demo-phase-diagrams.ipynb)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# here's an example of the phase diagrams that I create in pynamical-demo-phase-diagrams.ipynb\n",
    "IPdisplay.Image(url='images/3d-logistic-map-attractor-1.png', width=500)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## For more info:\n",
    "  - [Read the journal article](http://geoffboeing.com/publications/nonlinear-chaos-fractals-prediction/)\n",
    "  - [Pynamical documentation](http://pynamical.readthedocs.org)\n",
    "  - [Chaos Theory and the Logistic Map](http://geoffboeing.com/2015/03/chaos-theory-logistic-map/)\n",
    "  - [Visualizing Chaos and Randomness with Phase Diagrams](http://geoffboeing.com/2015/04/visualizing-chaos-and-randomness/)\n",
    "  - [Animated 3D Plots in Python](http://geoffboeing.com/2015/04/animated-3d-plots-python/)"
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python (pynamical)",
   "language": "python",
   "name": "pynamical"
  },
  "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.10.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}