{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## verhulst mandelbrot "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Most examples work across multiple plotting backends, this example is also available for:\n",
    "\n",
    "* [Bokeh - verhulst_mandelbrot ](../bokeh/verhulst_mandelbrot.ipynb)\n",
    "\n",
    "Example showing how bifurcation diagram for the logistic map relates to the Mandelbrot set according to a linear transformation. Inspired by [this illustration](https://en.wikipedia.org/wiki/Mandelbrot_set#/media/File:Verhulst-Mandelbrot-Bifurcation.jpg) on Wikipedia."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from itertools import islice\n",
    "\n",
    "import numpy as np\n",
    "import holoviews as hv\n",
    "from holoviews import opts\n",
    "\n",
    "hv.extension('matplotlib')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Defining Mandelbrot and Logistic Map"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Area of the complex plane\n",
    "bounds = (-2,-1.4,0.8,1.4)\n",
    "# Growth rates used in the logistic map\n",
    "growth_rates = np.linspace(0.9, 4, 1000)\n",
    "# Bifurcation points\n",
    "bifurcations = [1, 3, 3.4494, 3.5440, 3.5644, 3.7381, 3.7510, 3.8284, 3.8481]\n",
    "\n",
    "def mandelbrot_generator(h,w, maxit, bounds=bounds):\n",
    "    \"Generator that yields the mandlebrot set.\"\n",
    "    (l,b,r,t) = bounds\n",
    "    y,x = np.ogrid[b:t : h*1j, l:r:w*1j]\n",
    "    c = x+y*1j\n",
    "    z = c\n",
    "    divtime = maxit + np.zeros(z.shape, dtype=int)\n",
    "    for i in range(maxit):\n",
    "        z  = z**2 + c\n",
    "        diverge = z*np.conj(z) > 2**2\n",
    "        div_now = diverge & (divtime==maxit)\n",
    "        divtime[div_now] = i\n",
    "        z[diverge] = 2\n",
    "        yield divtime\n",
    "        \n",
    "def mandelbrot(h,w, n, maxit):\n",
    "    \"Returns the mandelbrot set computed to maxit\"\n",
    "    iterable =  mandelbrot_generator(h,w, maxit)\n",
    "    return next(islice(iterable, n, None))\n",
    "\n",
    "def mapping(r):\n",
    "    \"Linear mapping applied to the logistic bifurcation diagram\"\n",
    "    return (r /2.0) * ( 1 - (r/2.0))\n",
    "\n",
    "def logistic_map(gens=20, init=0.5, growth=0.5):\n",
    "    population = [init]\n",
    "    for gen in range(gens-1):\n",
    "        current = population[gen]\n",
    "        population.append(current * growth * (1 - current))\n",
    "    return population"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Plot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "bifurcation_diagram = hv.Points([(mapping(rate), pop) for rate in growth_rates for \n",
    "             (gen, pop) in enumerate(logistic_map(gens=200, growth=rate))\n",
    "             if gen>=100])  # Discard the first 100 generations to view attractors more easily\n",
    "\n",
    "img = hv.Image(mandelbrot(800,800, 45, 46).copy(), bounds=(-2, -1.4, 0.8, 1.4))\n",
    "\n",
    "vlines = hv.Overlay([hv.Curve([(mapping(pos),0), ((mapping(pos),1.4))]) for pos in bifurcations])\n",
    "\n",
    "(img * bifurcation_diagram * hv.HLine(0) * vlines).opts(\n",
    "    opts.Curve(color='teal', linewidth=1),\n",
    "    opts.HLine(color='k', linestyle='--'),\n",
    "    opts.Image(cmap='Reds', fig_size=200, logz=True, xaxis=None, yaxis=None),\n",
    "    opts.Points(s=1, color='g'))"
   ]
  }
 ],
 "metadata": {
  "language_info": {
   "name": "python",
   "pygments_lexer": "ipython3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}