samplecode/numpy/julia.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "33f94c77",
   "metadata": {},
   "source": [
    "# Filled Julia sets of quadratic polynomials\n",
    "\n",
    "A polished version of the notebook I started in Lecture 24."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 82,
   "id": "f019ac36",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import PIL.Image\n",
    "import time"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9c7eae98",
   "metadata": {},
   "source": [
    "## Reminders about complex numbers in python"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "94863ba2",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1j"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# In mathematics the square root of -1 is called i\n",
    "1j"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "f2ca164a",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(-1+0j)"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "1j*1j"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "fee71ad9",
   "metadata": {},
   "outputs": [],
   "source": [
    "z = 4-5j"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "77914495",
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "complex"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "type(z)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "d526d38d",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "6.4031242374328485"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "abs(z)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d443af98",
   "metadata": {},
   "source": [
    "## A few quadratic polynomials"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 93,
   "id": "54f8cb8b",
   "metadata": {},
   "outputs": [],
   "source": [
    "def k(z):\n",
    "    return z*z\n",
    "\n",
    "def q(z):\n",
    "    return z*z - 0.1\n",
    "\n",
    "def f(z):\n",
    "    return z*z - 1\n",
    "\n",
    "def g(z):\n",
    "    return z*z - 0.52 + 0.53j"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "864f1a78",
   "metadata": {},
   "outputs": [],
   "source": [
    "def orbit(func,a,n):\n",
    "    \"Compute the first `n` terms in the orbit of `a` under function `func`\"\n",
    "    L = [a]\n",
    "    for _ in range(n-1):\n",
    "        a = func(a)\n",
    "        L.append(a)\n",
    "    return L"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "a7e72b09",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[2,\n",
       " 3,\n",
       " 8,\n",
       " 63,\n",
       " 3968,\n",
       " 15745023,\n",
       " 247905749270528,\n",
       " 61457260521381894004129398783,\n",
       " 3776994870793005510047522464634252677140721938309041881088,\n",
       " 14265690253996672387291309349232388828298289458234016200317876247121873778287073518355813134107244701354409532063743]"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "orbit(f,2,10) # unbounded"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "920074e5",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[1.5,\n",
       " 1.25,\n",
       " 0.5625,\n",
       " -0.68359375,\n",
       " -0.5326995849609375,\n",
       " -0.7162311521824449,\n",
       " -0.4870129366434074,\n",
       " -0.7628183995419644,\n",
       " -0.418108089320236,\n",
       " -0.8251856256449815,\n",
       " -0.31906868322890036,\n",
       " -0.8981951753825757,\n",
       " -0.19324542691946411,\n",
       " -0.962656204974714,\n",
       " -0.07329303102368134,\n",
       " -0.9946281316033617,\n",
       " -0.010714879823205736,\n",
       " -0.9998851913503742,\n",
       " -0.00022960411822547933,\n",
       " -0.9999999472819489,\n",
       " -1.0543609951785271e-07,\n",
       " -0.9999999999999889,\n",
       " -2.220446049250313e-14,\n",
       " -1.0,\n",
       " 0.0,\n",
       " -1.0,\n",
       " 0.0,\n",
       " -1.0,\n",
       " 0.0,\n",
       " -1.0]"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "orbit(f,1.5,30) # bounded"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "id": "96ef14af",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[1.62,\n",
       " 1.6244000000000005,\n",
       " 1.6386753600000015,\n",
       " 1.6852569354711342,\n",
       " 1.8400909385535584,\n",
       " 2.3859346621469157,\n",
       " 4.692684212034117,\n",
       " 21.021285113874264,\n",
       " 440.8944278387917,\n",
       " 194386.89649929552,\n",
       " 37786265529.62783,\n",
       " 1.4278018626755405e+21,\n",
       " 2.0386181590597428e+42,\n",
       " 4.155963998448135e+84,\n",
       " 1.727203675639701e+169,\n",
       " inf,\n",
       " inf,\n",
       " inf,\n",
       " inf,\n",
       " inf]"
      ]
     },
     "execution_count": 32,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "orbit(f,1.62,20) # unbounded"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "id": "ec3ef3aa",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[(-0.5+0.5j),\n",
       " (-1-0.5j),\n",
       " (-0.25+1j),\n",
       " (-1.9375-0.5j),\n",
       " (2.50390625+1.9375j),\n",
       " (1.5156402587890625+9.70263671875j),\n",
       " (-92.84399390197359+29.41141365468502j),\n",
       " (7753.975950502717-5461.3462200079975j),\n",
       " (30297839.50617887-84694294.49462186j),\n",
       " (-6255164441199562-5132108283374204j),\n",
       " (1.2788546754169306e+31+6.42043624850961e+31j),\n",
       " (-3.958653234034041e+63+1.6421609829245704e+63j),\n",
       " (1.2974242733488377e+127-1.300149177171774e+127j),\n",
       " (-7.078137825681485e+251-3.3736902028743566e+254j),\n",
       " (nan+infj),\n",
       " (nan+nanj),\n",
       " (nan+nanj),\n",
       " (nan+nanj),\n",
       " (nan+nanj),\n",
       " (nan+nanj)]"
      ]
     },
     "execution_count": 35,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "orbit(f,-0.5+0.5j,20) # unbounded"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f7c52456",
   "metadata": {},
   "source": [
    "Using this function we can test individual points for bounded/unbounded orbits.  But how can we see all the answers at once?"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1c1538ab",
   "metadata": {},
   "source": [
    "## Goal\n",
    "\n",
    "Make a picture of (part of) the complex plane, where each point is colored black or white.\n",
    "* **white** for points where the orbit of `f` is unbounded\n",
    "* **black** for points where the orbit of `f` is bounded\n",
    "\n",
    "We can't actually tell an orbit is bounded with certainty, so we'll settle for this:\n",
    "* **white** for points where we could tell the orbit of `f` is unbounded, because some time in the first 1000 terms of that orbit we saw a value with absolute value greater than or equal to 2\n",
    "* **black** for points whose orbits are either bounded, or stay small for so long that we couldn't tell after 1000 terms\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7b23857b",
   "metadata": {},
   "source": [
    "## Grid of complex numbers"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3c9bad0f",
   "metadata": {},
   "source": [
    "As a first step, we need to make a matrix that stores a grid of complex numbers $x + iy$, evenly spaced in x and y coordinates."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "id": "c74eb84b",
   "metadata": {},
   "outputs": [],
   "source": [
    "size = 500 # How many horizontal and vertical grid points\n",
    "radius = 2.5 # We'll look at the square |x| <= radius and |y| <= radius\n",
    "xvals = np.linspace(-radius,radius,size) # row of all the real parts we need to use\n",
    "yvals = np.linspace(-radius,radius,size) # row of all the imaginary parts we need to use"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 52,
   "id": "514e3a72",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Combine a row of real parts and a column of imaginary parts into a matrix of complex numbers\n",
    "# [[ -2.5+2.5j, ..., 2.5+2.5j], ... [-2.5-2.5j, ..., 2.5-2.5j] ]\n",
    "# row            col\n",
    "zz = xvals + 1j * yvals[::-1].reshape((size,1))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8c2de5b2",
   "metadata": {},
   "source": [
    "The previous line is like taking\n",
    "$$\n",
    "(-2,-1,0,1,2) + \n",
    "\\begin{pmatrix}\n",
    "2i\\\\\n",
    "i\\\\\n",
    "0i\\\\\n",
    "-i\\\\\n",
    "-2i\n",
    "\\end{pmatrix}\n",
    "$$\n",
    "and having numpy broadcasating turn it into\n",
    "$$\n",
    "\\begin{pmatrix}\n",
    "-2 & -1 & 0 & 1 & 2\\\\\n",
    "-2 & -1 & 0 & 1 & 2\\\\\n",
    "-2 & -1 & 0 & 1 & 2\\\\\n",
    "-2 & -1 & 0 & 1 & 2\\\\\n",
    "-2 & -1 & 0 & 1 & 2\n",
    "\\end{pmatrix}\n",
    "+ \n",
    "\\begin{pmatrix}\n",
    "2i & 2i & 2i & 2i & 2i\\\\\n",
    "i & i & i & i & i\\\\\n",
    "0i & 0i & 0i & 0i & 0i\\\\\n",
    "-1i & -1i & -1i & -1i & -1i\\\\\n",
    "-2i & -2i & -2i & -2i & -2i\n",
    "\\end{pmatrix}\n",
    "$$\n",
    "which gives\n",
    "$$\n",
    "\\begin{pmatrix}\n",
    "-2+2i & -1+2i & 2i & 1+2i & 2+2i\\\\\n",
    "-2+i & -1+i & i & 1+i & 2+i\\\\\n",
    "-2 & -1 & 0 & 1 & 2\\\\\n",
    "-2-i & -1-i & 0-i & 1-i & 2-i\\\\\n",
    "-2-2i & -1-2i & -2i & 1-2i & 2-2i\n",
    "\\end{pmatrix}\n",
    "$$\n",
    "which is a matrix giving a rectangular grid of points in the complex plane."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "id": "d62be12b",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(500, 500)"
      ]
     },
     "execution_count": 53,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "## Confirm we have a (size)x(size) matrix now\n",
    "zz.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e96e738f",
   "metadata": {},
   "source": [
    "Let's classify the orbits of the points in zz"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "id": "9890857d",
   "metadata": {},
   "outputs": [],
   "source": [
    "def is_unbounded(a):\n",
    "    \"\"\"\n",
    "    Return True if the orbit of `a` under `f` can be seen to be unbounded\n",
    "    by examining the first 1000 points in it.  Otherwise, return False.\n",
    "    \"\"\"\n",
    "    for _ in range(1000):\n",
    "        if abs(a) >= 2:\n",
    "            return True\n",
    "        a = f(a)\n",
    "    \n",
    "    return False"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 55,
   "id": "16362920",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "True"
      ]
     },
     "execution_count": 55,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "is_unbounded(2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 56,
   "id": "2e1f62a1",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "True"
      ]
     },
     "execution_count": 56,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "is_unbounded(1.9)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 57,
   "id": "92ba75e8",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "False"
      ]
     },
     "execution_count": 57,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "is_unbounded(1.25)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "50aa8e19",
   "metadata": {},
   "source": [
    "Let's apply the `is_unbounded` function to every entry of the matrix "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f2912a98",
   "metadata": {},
   "source": [
    "Let's apply `is_unbounded` to every element of matrix `zz`, storing the results in a matrix of booleans."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 58,
   "id": "4cdaa143",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Where I'll store the orbit classification results\n",
    "imgdata = np.zeros_like(zz,dtype=\"uint8\")\n",
    "# black = 0, white = 255"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 59,
   "id": "f7eb3085",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(500, 500)"
      ]
     },
     "execution_count": 59,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "imgdata.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 85,
   "id": "99f32b84",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2.2318356037139893 seconds\n"
     ]
    }
   ],
   "source": [
    "t0 = time.time()\n",
    "\n",
    "for i in range(zz.shape[0]):\n",
    "    for j in range(zz.shape[1]):\n",
    "        if is_unbounded(zz[i,j]):\n",
    "            # orbit is not bounded, so put a white pixel at that point\n",
    "            imgdata[i,j] = 255\n",
    "\n",
    "            t1 = time.time()\n",
    "print(t1-t0,\"seconds\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 86,
   "id": "be9cdd48",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAfQAAAH0CAAAAADuvYBWAAAK8UlEQVR4nO2d23bjKhYAYdb5/19mHpw4ko1sXQDtpqoeOo5jIdhlLpKAziUJjf/dnQEZj9KBKB2I0oEoHYjSgSgdiNKBKB2I0oEoHYjSgSgdiNKBKB2I0oEoHYjSgSgdiNKBKB2I0oEoHYjSgSgdiNKBKB2I0oEoHYjSgSgdiNKBKB2I0oEoHYjSgSgdiNKBKB2I0oEoHYjSgSgdiNKBKB2I0oEoHYjSgSgdiNKBKB2I0oEoHYjSgSgdiNKBKB2I0oEoHYjSgSgdiNKBKB2I0oEoHYjSgSgdiNKBKB2I0oEoHYjSgSgdiNKBKB2I0oEoHYjSgSgdiNKBKB2I0oEoHYjSgSgdiNKBKB2I0oEoHYjSgSgdiNKBKB2I0oEoHYjSgSgdiNKBKB2I0oEoHYjSgSgdiNKBKB2I0oEoHYjSgSgdiNKBKB2I0oEoHYjSgSgdiNKBKB2I0oEoHYjSgSgdiNKBKB2I0oEoHYjSgSgdiNKBKB2I0oEoHYjSgSgdiNKBKB2I0oEoHYjSgSgdiNKBKB2I0oEoHYjSgSgdiNKBKB2I0oEoHYjSgSgdiNKBKB2I0oEoHYjSgSgdiNKBKB2I0oEoHYjSgSgdiNKBKB2I0oEoHYjSgSgdiNKBKB2I0oEoHQhXes535+A2sNJzSljrWOlksNJLSuXuPNwFVjoZpQPJ0DbuMYiDFh5a0/PqBw2m9Pz2AgVPes5L1Zl4kwbYp79K5kWAV9NfJfOcE6WL0oEoHQhQev74KwHg6H3tmVh+YE0XpPRSfckBKZ0OU3p5e4GCKf1XNtM5VfpDN9Q5Vrpz5IQF8uZMcroUjrycLgW8C5v+uzsDY8mppKXovHiXA6Sm/2ou72vYcv5r5iHVntGnP71uWS3rj03OjDX93ezXmpxfPvgpsX+fKWt6rVB75NUOmzE+cw7k3tvpXRX21fCMtTyliZv3vHpjp778cliaU/yE0teacn5768uxL8P7Ca3P2Get77Ydl7Y+csIATVjT15pOVNTVIRM6n3MgV361nWyaF+PAGZ1P2bw/uN4XTxuaGZv3lFKT8deEQ7gH4aRPGOlwRQonvRFNAh3OViMm7dNb6Zo0OtGKlVO6Hut2VbRJTqLFOFrzHq1ZjpafJsx3nd44wjM+Y49W06/TXNF0zieULl+JJv30ypPn/NbmHWj+S/s4IVfShBu91zvRHdn8Oa7HoKmknX17JZsRxwTRanpKpTpR7avLnht/7k27ls1ace4mnvTFvIcjmzmW55G9cnTE3qkiDCNk8/6k/GRw52qEngHemYGSHlkOva1NMOm5VMSVvbm8XXrKpZaJEi7Kd+dgxaa2Hdns3ZBeyUKsKEfq0y9ebnUO7KXkO1xJXiCS9E+Byd/D1jms30//adAWyXkI6fl37ejnT31ejdY/rF/OXBuNvH8qxGg+gvSfi5s9H6p+atSihIsZODQBvytBnrLtW3W0sVFMGbe+vHqmsnNvgxjGU5TR+/5wVLM7MpoXMxAi3DGa9/1Uw1uGhbJ+pjBVeCchavqxoJ1ch9yGyyePEO8INf2wsrcdRFrlZMepP+fkcAK3MFB6q+JGelhZmhZrFOOkb9xdOXGvat0lDb7Z9XK6b1fnXxP48nYXxkl/mloX7njUXpeVDm4wX79yFxP4S2Fc8zWweS/VHvCUsggdY0qpXe7z0K3sBo7e61/ps/42Go4RXD51PQAz1vRceXX+VnSIodzZTCxLXQ1LZ0bV9JcSXd+i8c5NHrvlfpCMm67Tr6u686HVvseCexK6nMIpOkvPLz9f/nCl0Dmn21r5ki5+3zaKvxmutvSu6XnrseP1cvWa5b6DNt+3rah0b7769uk/0ws2HpNcVjb6Er3lqb8EpauWjs/Tv41Lw1xs38JG6X8rYdd7zREeuMhgOtb0v4cR9bbwevP+L7PZvK9/dqHzdXreHm9dHofdOJDrl/ndayWv0Lt5L1trddsMfu+8ZGuSSu2t7hOBbrojl0r97QOUi8df4urJN4o/9x2566V7pDBudtzq3GWRhUsJXU7hFOOesi2+1aX25jEaJHGaXtkfpuKGSRTLwp2uqSEG/qcfES4jUHnVm4GLHcqHofzhtNbJDmWpqVVZxs77GzoxslavT5U1xNP0lNrlvgy9aXHDJIqXKWIn4raeeTL4G7BuWc7k/mWS4O+L+SZRfDrV0dK+bexw47Km6tYZRxJYpHQmO6cYOjGyWTohhnEpNe2L55wYuY3LmsYSQvqRuNXzO8r65bOHCPe/9mi1vmp0WE2vnymEyQNE2JRgv7LNzcXGXPFsnaW+kVj9oxG+IRGkp52Ds8dTjvegDeyiqt+6x+O+XfmPMQKN0byXtKON3HpIu/2HxlzMwK5CDiHGQO7B52rwddvF3rXo2+m/XLEHCnSMmv7Dh7iU77spd47q99N/enoUSHks6Rd3ye5c0y/OjopkPZL0lDa6zH0h6zuhYlfqG7U9kvCUYvXpD5Y1qiS3/u5APOk/97MfkSvLt74f1s962ZuF31fPIsRzHuU6fcXzCuhIuPLjoD7aHxk5UEEWRYhxab4iWp+eqjvu7OjUy+pHW/amffK/nxlOuOb97Gz2302iOzTxJR+u6X/cODt/k2jSG6zZbG29QWaCBTle8y7dmU968+Y9Xp98lYij92s0fpYVrGVuQrSa3iTGDUVFy08Tog3kGtGqrs8ZnWg1vZGuRrLaJBNuUBBOeiOmbJZbEU76hIEOV6Rw0lvRINLhZLVivku2lP7ufV7bvDXiHdQmTCl9uVbyjPeF6SmvbmZs3i8vK91MbBImlJ5ffjs2jeptpfiE1ieUvnZcDs3GKL9HbCU3BVP2WZUR2L76WjloxvhMOZCrmNo1t+L1sDDrkBozY02vLjd7LITbPKYsPvUtsX+eCfv06kS1rb+8/KE2PW8+ZpT+znPGYnWG5ee/TsiMrdcH8ns/XerN+szMOZDbpPz8+9fezzpC/wSspj8JOUt1FFTpvDZ9AWMgJyuw0vOkN172QJWen/8AgUrPqx80mNLz2wsUTOlwkNJz9SUHpHQ6wJszN/1naIEA1vTy8VcCQOmidCBKBwKUnj/+SoAn/eWRasjd/ToDvGRLy9qNLD6vpqe0UI10DpXedVPR+DCb95ScLiUslA4EK905coICK71gx+7k0TsYbE0no3QgSgeidCBKB6J0IEoHonQgSgeidCBKB6J0IEoHonQgSgeidCBKB6J0IEoHonQgSgeidCBKB6J0IEoHonQgSgeidCBKB6J0IEoHonQgSgeidCBKB6J0IEoHonQgSgeidCBKB6J0IEoHonQgSgeidCBKB6J0IEoHonQgSgeidCBKB6J0IEoHonQgSgeidCBKB6J0IEoHonQgSgeidCBKB6J0IEoHonQgSgeidCBKB6J0IEoHonQgSgeidCBKB6J0IEoHonQgSgeidCBKB6J0IEoHonQgSgeidCBKB6J0IEoHonQgSgeidCBKB6J0IEoHonQgSgeidCBKB6J0IEoHonQgSgeidCBKB6J0IEoHonQgSgeidCBKB6J0IEoHonQgSgeidCBKB6J0IEoHonQgSgeidCBKB6J0IEoHonQgSgeidCBKB6J0IEoHonQgSgeidCBKB6J0IEoHonQgSgeidCBKB6J0IEoHonQgSgeidCBKB6J0IEoHonQgSgeidCBKB6J0IEoHonQgSgeidCBKB6J0IEoHonQgSgeidCBKB6J0IEoHonQgSgeidCBKB/J/Vk3K+5Li/XcAAAAASUVORK5CYII=",
      "text/plain": [
       "<PIL.Image.Image image mode=L size=500x500>"
      ]
     },
     "execution_count": 86,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "PIL.Image.fromarray(imgdata)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "eecb345d",
   "metadata": {},
   "source": [
    "The center is black, because the orbit of 0 is bounded.  The black region stretches pretty far to the right, since for example we say the orbit of 1.5 is also unbounded.  But it doesn't go very high, since for example:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 63,
   "id": "2957bfbb",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "True"
      ]
     },
     "execution_count": 63,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "is_unbounded(1j)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "adefd0c1",
   "metadata": {},
   "source": [
    "the point $i$ or `1j` is already unbounded."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b3867233",
   "metadata": {},
   "source": [
    "## A more `numpy`-natural way to do this"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ee331887",
   "metadata": {},
   "source": [
    "`f` is agnostic about the type of its argument.  Since numpy arrays support squaring and addition, we can just apply `f` to the entire array!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 64,
   "id": "5adf3ac6",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[-1.       -12.5j       , -1.0499998-12.4498998j ,\n",
       "        -1.0997988-12.3997996j , ..., -1.0997988+12.3997996j ,\n",
       "        -1.0499998+12.4498998j , -1.       +12.5j       ],\n",
       "       [-0.9500002-12.4498998j , -1.       -12.4000004j ,\n",
       "        -1.049799 -12.350101j  , ..., -1.049799 +12.350101j  ,\n",
       "        -1.       +12.4000004j , -0.9500002+12.4498998j ],\n",
       "       [-0.9002012-12.3997996j , -0.950201 -12.350101j  ,\n",
       "        -1.       -12.30040241j, ..., -1.       +12.30040241j,\n",
       "        -0.950201 +12.350101j  , -0.9002012+12.3997996j ],\n",
       "       ...,\n",
       "       [-0.9002012+12.3997996j , -0.950201 +12.350101j  ,\n",
       "        -1.       +12.30040241j, ..., -1.       -12.30040241j,\n",
       "        -0.950201 -12.350101j  , -0.9002012-12.3997996j ],\n",
       "       [-0.9500002+12.4498998j , -1.       +12.4000004j ,\n",
       "        -1.049799 +12.350101j  , ..., -1.049799 -12.350101j  ,\n",
       "        -1.       -12.4000004j , -0.9500002-12.4498998j ],\n",
       "       [-1.       +12.5j       , -1.0499998+12.4498998j ,\n",
       "        -1.0997988+12.3997996j , ..., -1.0997988-12.3997996j ,\n",
       "        -1.0499998-12.4498998j , -1.       -12.5j       ]])"
      ]
     },
     "execution_count": 64,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "f(zz)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7291218c",
   "metadata": {},
   "source": [
    "Idea: Let's make a copy of the grid `zz` and keep replacing that copy with `f` applied to it.  Points with unbounded orbits will overflow, because the numbers get so big so fast.  Any points that don't overflow are bounded (probably)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 66,
   "id": "72d1ab05",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Users\\ddumas\\AppData\\Local\\Temp\\ipykernel_7760\\3119929435.py:5: RuntimeWarning: overflow encountered in multiply\n",
      "  return z*z - 1\n",
      "C:\\Users\\ddumas\\AppData\\Local\\Temp\\ipykernel_7760\\3119929435.py:5: RuntimeWarning: invalid value encountered in multiply\n",
      "  return z*z - 1\n"
     ]
    }
   ],
   "source": [
    "iterates = np.copy(zz)\n",
    "for _ in range(1000):\n",
    "    iterates = f(iterates)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 67,
   "id": "fc7b63d4",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(nan+nanj)"
      ]
     },
     "execution_count": 67,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "iterates[0,0] # unbounded point -> overflow -> \"NAN\" value (not a number)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 71,
   "id": "d3e133ab",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0j"
      ]
     },
     "execution_count": 71,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "iterates[251,253] # bounded point -> value that is not NAN"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f351ad43",
   "metadata": {},
   "source": [
    "`numpy` has a ufunc `isnan` that tests whether a value is NAN.  In this case, that's equivalent to us thinking it is a point with unbounded `f`-orbit."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 72,
   "id": "f7dce7b3",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ True,  True,  True, ...,  True,  True,  True],\n",
       "       [ True,  True,  True, ...,  True,  True,  True],\n",
       "       [ True,  True,  True, ...,  True,  True,  True],\n",
       "       ...,\n",
       "       [ True,  True,  True, ...,  True,  True,  True],\n",
       "       [ True,  True,  True, ...,  True,  True,  True],\n",
       "       [ True,  True,  True, ...,  True,  True,  True]])"
      ]
     },
     "execution_count": 72,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.isnan(iterates) # is_unbounded for each point in the grid"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a69ed7df",
   "metadata": {},
   "source": [
    "So now we can put that together with the image conversion step into a 4-line cell that computes the filled julia set!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 84,
   "id": "850a6940",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Users\\ddumas\\AppData\\Local\\Temp\\ipykernel_7760\\3119929435.py:5: RuntimeWarning: overflow encountered in multiply\n",
      "  return z*z - 1\n",
      "C:\\Users\\ddumas\\AppData\\Local\\Temp\\ipykernel_7760\\3119929435.py:5: RuntimeWarning: invalid value encountered in multiply\n",
      "  return z*z - 1\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2.1183784008026123 seconds\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAfQAAAH0CAAAAADuvYBWAAAK8UlEQVR4nO2d23bjKhYAYdb5/19mHpw4ko1sXQDtpqoeOo5jIdhlLpKAziUJjf/dnQEZj9KBKB2I0oEoHYjSgSgdiNKBKB2I0oEoHYjSgSgdiNKBKB2I0oEoHYjSgSgdiNKBKB2I0oEoHYjSgSgdiNKBKB2I0oEoHYjSgSgdiNKBKB2I0oEoHYjSgSgdiNKBKB2I0oEoHYjSgSgdiNKBKB2I0oEoHYjSgSgdiNKBKB2I0oEoHYjSgSgdiNKBKB2I0oEoHYjSgSgdiNKBKB2I0oEoHYjSgSgdiNKBKB2I0oEoHYjSgSgdiNKBKB2I0oEoHYjSgSgdiNKBKB2I0oEoHYjSgSgdiNKBKB2I0oEoHYjSgSgdiNKBKB2I0oEoHYjSgSgdiNKBKB2I0oEoHYjSgSgdiNKBKB2I0oEoHYjSgSgdiNKBKB2I0oEoHYjSgSgdiNKBKB2I0oEoHYjSgSgdiNKBKB2I0oEoHYjSgSgdiNKBKB2I0oEoHYjSgSgdiNKBKB2I0oEoHYjSgSgdiNKBKB2I0oEoHYjSgSgdiNKBKB2I0oEoHYjSgSgdiNKBKB2I0oEoHYjSgSgdiNKBKB2I0oEoHYjSgSgdiNKBKB2I0oEoHYjSgSgdiNKBKB2I0oEoHQhXes535+A2sNJzSljrWOlksNJLSuXuPNwFVjoZpQPJ0DbuMYiDFh5a0/PqBw2m9Pz2AgVPes5L1Zl4kwbYp79K5kWAV9NfJfOcE6WL0oEoHQhQev74KwHg6H3tmVh+YE0XpPRSfckBKZ0OU3p5e4GCKf1XNtM5VfpDN9Q5Vrpz5IQF8uZMcroUjrycLgW8C5v+uzsDY8mppKXovHiXA6Sm/2ou72vYcv5r5iHVntGnP71uWS3rj03OjDX93ezXmpxfPvgpsX+fKWt6rVB75NUOmzE+cw7k3tvpXRX21fCMtTyliZv3vHpjp778cliaU/yE0teacn5768uxL8P7Ca3P2Get77Ydl7Y+csIATVjT15pOVNTVIRM6n3MgV361nWyaF+PAGZ1P2bw/uN4XTxuaGZv3lFKT8deEQ7gH4aRPGOlwRQonvRFNAh3OViMm7dNb6Zo0OtGKlVO6Hut2VbRJTqLFOFrzHq1ZjpafJsx3nd44wjM+Y49W06/TXNF0zieULl+JJv30ypPn/NbmHWj+S/s4IVfShBu91zvRHdn8Oa7HoKmknX17JZsRxwTRanpKpTpR7avLnht/7k27ls1ace4mnvTFvIcjmzmW55G9cnTE3qkiDCNk8/6k/GRw52qEngHemYGSHlkOva1NMOm5VMSVvbm8XXrKpZaJEi7Kd+dgxaa2Hdns3ZBeyUKsKEfq0y9ebnUO7KXkO1xJXiCS9E+Byd/D1jms30//adAWyXkI6fl37ejnT31ejdY/rF/OXBuNvH8qxGg+gvSfi5s9H6p+atSihIsZODQBvytBnrLtW3W0sVFMGbe+vHqmsnNvgxjGU5TR+/5wVLM7MpoXMxAi3DGa9/1Uw1uGhbJ+pjBVeCchavqxoJ1ch9yGyyePEO8INf2wsrcdRFrlZMepP+fkcAK3MFB6q+JGelhZmhZrFOOkb9xdOXGvat0lDb7Z9XK6b1fnXxP48nYXxkl/mloX7njUXpeVDm4wX79yFxP4S2Fc8zWweS/VHvCUsggdY0qpXe7z0K3sBo7e61/ps/42Go4RXD51PQAz1vRceXX+VnSIodzZTCxLXQ1LZ0bV9JcSXd+i8c5NHrvlfpCMm67Tr6u686HVvseCexK6nMIpOkvPLz9f/nCl0Dmn21r5ki5+3zaKvxmutvSu6XnrseP1cvWa5b6DNt+3rah0b7769uk/0ws2HpNcVjb6Er3lqb8EpauWjs/Tv41Lw1xs38JG6X8rYdd7zREeuMhgOtb0v4cR9bbwevP+L7PZvK9/dqHzdXreHm9dHofdOJDrl/ndayWv0Lt5L1trddsMfu+8ZGuSSu2t7hOBbrojl0r97QOUi8df4urJN4o/9x2566V7pDBudtzq3GWRhUsJXU7hFOOesi2+1aX25jEaJHGaXtkfpuKGSRTLwp2uqSEG/qcfES4jUHnVm4GLHcqHofzhtNbJDmWpqVVZxs77GzoxslavT5U1xNP0lNrlvgy9aXHDJIqXKWIn4raeeTL4G7BuWc7k/mWS4O+L+SZRfDrV0dK+bexw47Km6tYZRxJYpHQmO6cYOjGyWTohhnEpNe2L55wYuY3LmsYSQvqRuNXzO8r65bOHCPe/9mi1vmp0WE2vnymEyQNE2JRgv7LNzcXGXPFsnaW+kVj9oxG+IRGkp52Ds8dTjvegDeyiqt+6x+O+XfmPMQKN0byXtKON3HpIu/2HxlzMwK5CDiHGQO7B52rwddvF3rXo2+m/XLEHCnSMmv7Dh7iU77spd47q99N/enoUSHks6Rd3ye5c0y/OjopkPZL0lDa6zH0h6zuhYlfqG7U9kvCUYvXpD5Y1qiS3/u5APOk/97MfkSvLt74f1s962ZuF31fPIsRzHuU6fcXzCuhIuPLjoD7aHxk5UEEWRYhxab4iWp+eqjvu7OjUy+pHW/amffK/nxlOuOb97Gz2302iOzTxJR+u6X/cODt/k2jSG6zZbG29QWaCBTle8y7dmU968+Y9Xp98lYij92s0fpYVrGVuQrSa3iTGDUVFy08Tog3kGtGqrs8ZnWg1vZGuRrLaJBNuUBBOeiOmbJZbEU76hIEOV6Rw0lvRINLhZLVivku2lP7ufV7bvDXiHdQmTCl9uVbyjPeF6SmvbmZs3i8vK91MbBImlJ5ffjs2jeptpfiE1ieUvnZcDs3GKL9HbCU3BVP2WZUR2L76WjloxvhMOZCrmNo1t+L1sDDrkBozY02vLjd7LITbPKYsPvUtsX+eCfv06kS1rb+8/KE2PW8+ZpT+znPGYnWG5ee/TsiMrdcH8ns/XerN+szMOZDbpPz8+9fezzpC/wSspj8JOUt1FFTpvDZ9AWMgJyuw0vOkN172QJWen/8AgUrPqx80mNLz2wsUTOlwkNJz9SUHpHQ6wJszN/1naIEA1vTy8VcCQOmidCBKBwKUnj/+SoAn/eWRasjd/ToDvGRLy9qNLD6vpqe0UI10DpXedVPR+DCb95ScLiUslA4EK905coICK71gx+7k0TsYbE0no3QgSgeidCBKB6J0IEoHonQgSgeidCBKB6J0IEoHonQgSgeidCBKB6J0IEoHonQgSgeidCBKB6J0IEoHonQgSgeidCBKB6J0IEoHonQgSgeidCBKB6J0IEoHonQgSgeidCBKB6J0IEoHonQgSgeidCBKB6J0IEoHonQgSgeidCBKB6J0IEoHonQgSgeidCBKB6J0IEoHonQgSgeidCBKB6J0IEoHonQgSgeidCBKB6J0IEoHonQgSgeidCBKB6J0IEoHonQgSgeidCBKB6J0IEoHonQgSgeidCBKB6J0IEoHonQgSgeidCBKB6J0IEoHonQgSgeidCBKB6J0IEoHonQgSgeidCBKB6J0IEoHonQgSgeidCBKB6J0IEoHonQgSgeidCBKB6J0IEoHonQgSgeidCBKB6J0IEoHonQgSgeidCBKB6J0IEoHonQgSgeidCBKB6J0IEoHonQgSgeidCBKB6J0IEoHonQgSgeidCBKB6J0IEoHonQgSgeidCBKB6J0IEoHonQgSgeidCBKB6J0IEoHonQgSgeidCBKB6J0IEoHonQgSgeidCBKB/J/Vk3K+5Li/XcAAAAASUVORK5CYII=",
      "text/plain": [
       "<PIL.Image.Image image mode=L size=500x500>"
      ]
     },
     "execution_count": 84,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "t0 = time.time()\n",
    "iterates = np.copy(zz)\n",
    "for _ in range(1000):\n",
    "    iterates = f(iterates)\n",
    "t1 = time.time()\n",
    "print(t1-t0,\"seconds\")\n",
    "PIL.Image.fromarray(255*np.isnan(iterates).astype(\"uint8\"))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "acdb7646",
   "metadata": {},
   "source": [
    "## An even better way"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c14b00e4",
   "metadata": {},
   "source": [
    "If a point has unbounded orbit, we sometimes know that after just a few terms.  But the code above still applies `f` 1000 times.  A better way would be to keep track of which points we still need to work on, and only iterate on them.  For the others, we don't do any calculations.\n",
    "\n",
    "We'll use a mask for this: An array of booleans `active` that tells whether each pixel is still supposed to be iterated.\n",
    "At each step, pixels where the value has gone outside $|z|<2$ will be inactivated by setting `active` to `False`.  But we'll only check the absolute value for active points; anything else is already inactive."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 83,
   "id": "d1d65d4a",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.24405121803283691 seconds\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAfQAAAH0CAAAAADuvYBWAAAK8UlEQVR4nO2d23bjKhYAYdb5/19mHpw4ko1sXQDtpqoeOo5jIdhlLpKAziUJjf/dnQEZj9KBKB2I0oEoHYjSgSgdiNKBKB2I0oEoHYjSgSgdiNKBKB2I0oEoHYjSgSgdiNKBKB2I0oEoHYjSgSgdiNKBKB2I0oEoHYjSgSgdiNKBKB2I0oEoHYjSgSgdiNKBKB2I0oEoHYjSgSgdiNKBKB2I0oEoHYjSgSgdiNKBKB2I0oEoHYjSgSgdiNKBKB2I0oEoHYjSgSgdiNKBKB2I0oEoHYjSgSgdiNKBKB2I0oEoHYjSgSgdiNKBKB2I0oEoHYjSgSgdiNKBKB2I0oEoHYjSgSgdiNKBKB2I0oEoHYjSgSgdiNKBKB2I0oEoHYjSgSgdiNKBKB2I0oEoHYjSgSgdiNKBKB2I0oEoHYjSgSgdiNKBKB2I0oEoHYjSgSgdiNKBKB2I0oEoHYjSgSgdiNKBKB2I0oEoHYjSgSgdiNKBKB2I0oEoHYjSgSgdiNKBKB2I0oEoHYjSgSgdiNKBKB2I0oEoHYjSgSgdiNKBKB2I0oEoHYjSgSgdiNKBKB2I0oEoHYjSgSgdiNKBKB2I0oEoHYjSgSgdiNKBKB2I0oEoHYjSgSgdiNKBKB2I0oEoHQhXes535+A2sNJzSljrWOlksNJLSuXuPNwFVjoZpQPJ0DbuMYiDFh5a0/PqBw2m9Pz2AgVPes5L1Zl4kwbYp79K5kWAV9NfJfOcE6WL0oEoHQhQev74KwHg6H3tmVh+YE0XpPRSfckBKZ0OU3p5e4GCKf1XNtM5VfpDN9Q5Vrpz5IQF8uZMcroUjrycLgW8C5v+uzsDY8mppKXovHiXA6Sm/2ou72vYcv5r5iHVntGnP71uWS3rj03OjDX93ezXmpxfPvgpsX+fKWt6rVB75NUOmzE+cw7k3tvpXRX21fCMtTyliZv3vHpjp778cliaU/yE0teacn5768uxL8P7Ca3P2Get77Ydl7Y+csIATVjT15pOVNTVIRM6n3MgV361nWyaF+PAGZ1P2bw/uN4XTxuaGZv3lFKT8deEQ7gH4aRPGOlwRQonvRFNAh3OViMm7dNb6Zo0OtGKlVO6Hut2VbRJTqLFOFrzHq1ZjpafJsx3nd44wjM+Y49W06/TXNF0zieULl+JJv30ypPn/NbmHWj+S/s4IVfShBu91zvRHdn8Oa7HoKmknX17JZsRxwTRanpKpTpR7avLnht/7k27ls1ace4mnvTFvIcjmzmW55G9cnTE3qkiDCNk8/6k/GRw52qEngHemYGSHlkOva1NMOm5VMSVvbm8XXrKpZaJEi7Kd+dgxaa2Hdns3ZBeyUKsKEfq0y9ebnUO7KXkO1xJXiCS9E+Byd/D1jms30//adAWyXkI6fl37ejnT31ejdY/rF/OXBuNvH8qxGg+gvSfi5s9H6p+atSihIsZODQBvytBnrLtW3W0sVFMGbe+vHqmsnNvgxjGU5TR+/5wVLM7MpoXMxAi3DGa9/1Uw1uGhbJ+pjBVeCchavqxoJ1ch9yGyyePEO8INf2wsrcdRFrlZMepP+fkcAK3MFB6q+JGelhZmhZrFOOkb9xdOXGvat0lDb7Z9XK6b1fnXxP48nYXxkl/mloX7njUXpeVDm4wX79yFxP4S2Fc8zWweS/VHvCUsggdY0qpXe7z0K3sBo7e61/ps/42Go4RXD51PQAz1vRceXX+VnSIodzZTCxLXQ1LZ0bV9JcSXd+i8c5NHrvlfpCMm67Tr6u686HVvseCexK6nMIpOkvPLz9f/nCl0Dmn21r5ki5+3zaKvxmutvSu6XnrseP1cvWa5b6DNt+3rah0b7769uk/0ws2HpNcVjb6Er3lqb8EpauWjs/Tv41Lw1xs38JG6X8rYdd7zREeuMhgOtb0v4cR9bbwevP+L7PZvK9/dqHzdXreHm9dHofdOJDrl/ndayWv0Lt5L1trddsMfu+8ZGuSSu2t7hOBbrojl0r97QOUi8df4urJN4o/9x2566V7pDBudtzq3GWRhUsJXU7hFOOesi2+1aX25jEaJHGaXtkfpuKGSRTLwp2uqSEG/qcfES4jUHnVm4GLHcqHofzhtNbJDmWpqVVZxs77GzoxslavT5U1xNP0lNrlvgy9aXHDJIqXKWIn4raeeTL4G7BuWc7k/mWS4O+L+SZRfDrV0dK+bexw47Km6tYZRxJYpHQmO6cYOjGyWTohhnEpNe2L55wYuY3LmsYSQvqRuNXzO8r65bOHCPe/9mi1vmp0WE2vnymEyQNE2JRgv7LNzcXGXPFsnaW+kVj9oxG+IRGkp52Ds8dTjvegDeyiqt+6x+O+XfmPMQKN0byXtKON3HpIu/2HxlzMwK5CDiHGQO7B52rwddvF3rXo2+m/XLEHCnSMmv7Dh7iU77spd47q99N/enoUSHks6Rd3ye5c0y/OjopkPZL0lDa6zH0h6zuhYlfqG7U9kvCUYvXpD5Y1qiS3/u5APOk/97MfkSvLt74f1s962ZuF31fPIsRzHuU6fcXzCuhIuPLjoD7aHxk5UEEWRYhxab4iWp+eqjvu7OjUy+pHW/amffK/nxlOuOb97Gz2302iOzTxJR+u6X/cODt/k2jSG6zZbG29QWaCBTle8y7dmU968+Y9Xp98lYij92s0fpYVrGVuQrSa3iTGDUVFy08Tog3kGtGqrs8ZnWg1vZGuRrLaJBNuUBBOeiOmbJZbEU76hIEOV6Rw0lvRINLhZLVivku2lP7ufV7bvDXiHdQmTCl9uVbyjPeF6SmvbmZs3i8vK91MbBImlJ5ffjs2jeptpfiE1ieUvnZcDs3GKL9HbCU3BVP2WZUR2L76WjloxvhMOZCrmNo1t+L1sDDrkBozY02vLjd7LITbPKYsPvUtsX+eCfv06kS1rb+8/KE2PW8+ZpT+znPGYnWG5ee/TsiMrdcH8ns/XerN+szMOZDbpPz8+9fezzpC/wSspj8JOUt1FFTpvDZ9AWMgJyuw0vOkN172QJWen/8AgUrPqx80mNLz2wsUTOlwkNJz9SUHpHQ6wJszN/1naIEA1vTy8VcCQOmidCBKBwKUnj/+SoAn/eWRasjd/ToDvGRLy9qNLD6vpqe0UI10DpXedVPR+DCb95ScLiUslA4EK905coICK71gx+7k0TsYbE0no3QgSgeidCBKB6J0IEoHonQgSgeidCBKB6J0IEoHonQgSgeidCBKB6J0IEoHonQgSgeidCBKB6J0IEoHonQgSgeidCBKB6J0IEoHonQgSgeidCBKB6J0IEoHonQgSgeidCBKB6J0IEoHonQgSgeidCBKB6J0IEoHonQgSgeidCBKB6J0IEoHonQgSgeidCBKB6J0IEoHonQgSgeidCBKB6J0IEoHonQgSgeidCBKB6J0IEoHonQgSgeidCBKB6J0IEoHonQgSgeidCBKB6J0IEoHonQgSgeidCBKB6J0IEoHonQgSgeidCBKB6J0IEoHonQgSgeidCBKB6J0IEoHonQgSgeidCBKB6J0IEoHonQgSgeidCBKB6J0IEoHonQgSgeidCBKB6J0IEoHonQgSgeidCBKB6J0IEoHonQgSgeidCBKB6J0IEoHonQgSgeidCBKB6J0IEoHonQgSgeidCBKB6J0IEoHonQgSgeidCBKB6J0IEoHonQgSgeidCBKB6J0IEoHonQgSgeidCBKB6J0IEoHonQgSgeidCBKB6J0IEoHonQgSgeidCBKB/J/Vk3K+5Li/XcAAAAASUVORK5CYII=",
      "text/plain": [
       "<PIL.Image.Image image mode=L size=500x500>"
      ]
     },
     "execution_count": 83,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "t0 = time.time()\n",
    "iterates = np.copy(zz)\n",
    "active = np.ones_like(iterates,dtype=\"bool\")  # Grid of True\n",
    "for _ in range(1000):\n",
    "    iterates[active] = f(iterates[active])  # Apply f to each pixel where we are supposed to proceed\n",
    "    active[active] &= np.abs(iterates[active]) < 2  # In future iterations, only continue on the pixels where abs < 2\n",
    "                                                     # (but only check this for points that remain active)\n",
    "t1 = time.time()\n",
    "print(t1-t0,\"seconds\")\n",
    "PIL.Image.fromarray(255*(1 - active.astype(\"uint8\")))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "401cb206",
   "metadata": {},
   "source": [
    "Note that mask behavior in numpy means `active[active]` is a way of saying \"the entries in array `active` that are `True` right now."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "aa175538",
   "metadata": {},
   "source": [
    "## Summary\n",
    "\n",
    "By using numpy conventions and features, we got the calculation time down by a factor of 10!"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6c19f904",
   "metadata": {},
   "source": [
    "## Some other nice Julia sets"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d205e6d2",
   "metadata": {},
   "source": [
    "### $k(z) = z^2$ (the unit disk)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 90,
   "id": "cd802f6b",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.538597583770752 seconds\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAfQAAAH0CAAAAADuvYBWAAAHuklEQVR4nO3Zy1IcSRAF0WRs/v+XNQsGPVpAd0Nl1b3h7jutFBnHokCmlx/LaP1z9QB2fqIDEx2Y6MBEByY6MNGBiQ5MdGCiAxMdmOjARAcmOjDRgYkOTHRgogMTHZjowEQHJjow0YGJDkx0YKIDEx2Y6MBEByY6MNGBiQ5MdGCiAxMdmOjARAcmOjDRgYkOTHRgogMTHZjowEQHJjow0YGJDkx0YKIDEx2Y6MBEByY6MNGBiQ5MdGCiAxMdmOjARAcmOjDRgYkOTHRgogMTHZjowEQHJjow0YGJDkx0YKIDEx2Y6MBEByY6MNGBiQ5MdGCiAxMdmOjARAcmOjDRgYkOTHRgogMTHZjowEQHJjow0YGJDkx0YKIDEx2Y6MBEByY6MNGBiQ5MdGCiAxMdmOjARAcmOjDRgYkOTHRgogMTHZjowEQHJjow0YGJDkx0YKIDEx2Y6MBEByY6MNGBiQ5MdGCiAxMdmOjARAcmOjDRgYkOTHRgogMTHZjowEQHJjow0YGJDkx0YKIDEx2Y6MBEB/bv1QOc38vNn39cMsWVvXCefIt9G2gTiKfe8/4VYx3jX/k4+FvzVzL6hc+DvzV7LXNf93Xx1wZvZujTviv+2tTljHzXMeRrDWUfiH6c+GsDNzTtSUeTrzWPfRj6DvK1prGPQt9FvtYs9kHoO8nXmsQ+Bn03+Vpz2Kf81+oZ5uf8JSc049LP05ixrgGvOPcAByxswOf95I/ugG98/aVfYdC+s/ZLv+Tu2o+9HP2i9ZerV3/er9x9896aL/3Se2s+9mL0i9derN6LfvnSLx/gy9WiB6w8YISvVfqLXMq+O7fXeekp5jmDPFUletCqg0Z5vEb0qEVHDfNghehhaw4b55H60OOWHDfQ3erQA1ccONLn1aHb92tDj7yqyKE+qQw9dL2hY31UF3rscmMHe7cq9ODVBo/2d03o0YuNHu6mJnQ7qCL08FsKH+/3etDjlxo/4M9q0AtWWjDiazXodlwt6BVXVDHkqkEvWWfJmCXodmQd6CUX1DJoB7odWgV6x/2stUpGbUCvWORbDcM2oNvBFaA33M5vFYxbgG5Hl49ecDl/lj9wProdXjx6/t38VfzI8eh2fOno8VfzXulDp6PbhsLR02/mg8LHDke3HYkOLBs9/DP5cdmDZ6PblqLRs+/l06JHj0a3PYkOLBk9+hN5r+Thk9FtU6IDC0ZP/kA+UPD4wei2K9GBiQ4sFz34Z+Jj5T4gF922JTow0YHFouf+RHy42CfEotu+RAcmOjDRgaWix/4S9Eypj0hFt42JDkx0YKIDEx2Y6MBEBxaKnvov3CcLfUYouu1MdGCiAxMdmOjARAcmOjDRgYkOTHRgogMTHZjowEQHJjqwUPQfVw9wTKHPCEW3nYkOTHRgogMTHZjowEQHlooe+i/c50p9RCq6bUx0YKIDEx1YLHrqL0FPFPuEWHTbl+jARAeWix77E/HRch+Qi27bEh2Y6MCC0XN/Jj5U8PjB6LYr0YElowd/IO+XPHwyum1KdGDR6MmfyDtFjx6NbnvKRo++l8/KHjwb3bYkOrBw9OzP5IeFjx2ObjtKRw+/mfdLHzod3TYUj55+Ne8UP3I8uh1fPnr83dyWP3A+uh1eAXr+5fxRwbgF6HZ0DegFt/OrhmEb0CsW+X8Vo1ag27F1oFfcz1otg3ag26GVoHdcUMuYJegd66wYcvWg24HVoBdcUcGIr9Wg5680fsCf9aCnLzV8vN8rQrejakKPvqXo4W5qQk9ebPBof1eFnrva2MHerQs9dbmhY31UGXrmeiOH+qQ2dDugOvTAqwoc6fPq0PNWHDfQ3frQ05YcNs4jFaJnrTlqmAdrRE9adNAoj1eJnrPqmEGe6qVz7LVerh5gtZK3XvqKWHjACF+rFv36lV8+wJfrRb966b3mzejXrr3YvPcXudeu+nWue2vNl74uW363eTv6NesvN2//vK91/ie+f2Ptl75OR+g3n3Dp68xjn7GuEa84i33IsgZ83tdaJ3EMMR9z6Wv/sQ/a1Jyn7GUftadJj9nHPmxLs56zh33cjqY96Hj2gRua96R1pPvM9Yx81VHsU5cz9F3r++6DNzP3aes77rPXMvp1a30Ffv5Kxr9wrWfgGetAvHKtdV8etAnOU9+6xQdugPdkm/Jfq/ZEogMTHZjowEQHJjow0YGJDkx0YKIDEx2Y6MBEByY6MNGBiQ5MdGCiAxMdmOjARAcmOjDRgYkOTHRgogMTHZjowEQHJjow0YGJDkx0YKIDEx2Y6MBEByY6MNGBiQ5MdGCiAxMdmOjARAcmOjDRgYkOTHRgogMTHZjowEQHJjow0YGJDkx0YKIDEx2Y6MBEByY6MNGBiQ5MdGCiAxMdmOjARAcmOjDRgYkOTHRgogMTHZjowEQHJjow0YGJDkx0YKIDEx2Y6MBEByY6MNGBiQ5MdGCiAxMdmOjARAcmOjDRgYkOTHRgogMTHZjowEQHJjow0YGJDkx0YKIDEx2Y6MBEByY6MNGBiQ5MdGCiAxMdmOjARAcmOjDRgYkOTHRgogMTHZjowEQHJjow0YGJDkx0YKIDEx2Y6MBEByY6MNGBiQ5MdGCiAxMdmOjARAcmOjDRgYkOTHRgogMTHZjowEQH9h+4gZ6fvM346QAAAABJRU5ErkJggg==",
      "text/plain": [
       "<PIL.Image.Image image mode=L size=500x500>"
      ]
     },
     "execution_count": 90,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "t0 = time.time()\n",
    "iterates = np.copy(zz)\n",
    "active = np.ones_like(iterates,dtype=\"bool\")  # Grid of True\n",
    "for _ in range(1000):\n",
    "    iterates[active] = k(iterates[active])  # Apply f to each pixel where we are supposed to proceed\n",
    "    active[active] &= np.abs(iterates[active]) < 2  # In future iterations, only continue on the pixels where abs < 2\n",
    "                                                     # (but only check this for points that remain active)\n",
    "t1 = time.time()\n",
    "print(t1-t0,\"seconds\")\n",
    "PIL.Image.fromarray(255*(1 - active.astype(\"uint8\")))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ea159cdf",
   "metadata": {},
   "source": [
    "### $q(z) = z^2 - 0.1$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 94,
   "id": "61ccf3b8",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.6459994316101074 seconds\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAfQAAAH0CAAAAADuvYBWAAAILElEQVR4nO3Yy24cRxAF0aKh//9leUFRlMQhOY+u6Xsz4mwMeGFkZSB7CL/8XKL57+wB9HxGBzI6kNGBjA5kdCCjAxkdyOhARgcyOpDRgYwOZHQgowMZHcjoQEYHMjqQ0YGMDmR0IKMDGR3I6EBGBzI6kNGBjA5kdCCjAxkdyOhARgcyOpDRgYwOZHQgowMZHcjoQEYHMjqQ0YGMDmR0IKMDGR3I6EBGBzI6kNGBjA5kdCCjAxkdyOhARgcyOpDRgYwOZHQgowMZHcjoQEYHMjqQ0YGMDmR0IKMDGR3I6EBGBzI6kNGBjA5kdCCjAxkdyOhARgcyOpDRgYwOZHQgowMZHcjoQEYHMjqQ0YGMDmR0IKMDGR3I6EBGBzI6kNGBjA5kdCCjAxkdyOhARgcyOpDRgYwOZHQgowMZHcjoQEYHMjqQ0YGMDmR0IKMDGR3I6EBGBzI6kNGBjA5kdCCjAxkdyOhARgcyOpDRgYwOZHQgowMZHcjoQEYHMjqQ0YGMDmR0IKMDGR3I6EBGBzI6kNGBjA5kdKAfZw9wjpdf//x56hRnecE9++XDv+GtAPbij8nXwmUHRb/c+x1oE5Cnflf8FWUZiHdel3wtSHZA9OuLvwJsZPYTbw3+ZvhWJj/v3uRrzc4+Nvojwd+M3c3Ehx0R/M3E/Uz8f+9HNj/2PxZiYPSDMw2sPu7zvqPRtB1Nu/Qtdznt2Gdd+r46s9Y06DV7D3LQogZ93jd/hAd948dc+jOajNnVjIc86wyHbGvCM5755R2xr/pHPP+3tn5l9X/InfD3Vf2fdO3RTwnQXr37837e9rvXVjz9uQdXvLjiz/vJH9nib3zrpWesvHR5pZee0TxljFtVXnrSrhv313jpSc2zhrlS36XnbbluhXWXntc8caSvtUWPXHDkUF8oix663tCxPtMVPXa5sYNd1BVdh6iKHnxPwaN91BQ9erHRw/2jKHr4WsPH+1NRdB2lJ3r8JcUP+FtPdB2mJnrBHRWM+KolesVCK4ZcNdFL1lkyZkf0kmW2DFoRvWOVa62SUSui61gN0Suu503DsA3RdbCC6A2384eCcfOjFyzxb/kD50fX4eKj59/NB/Ejp0ePX+Al6UOHR09f3yfCxw6Prh2yo4dfzOeyB8+Ori2io2ffy5eiR0+OHr247yQPHxw9eW1XCB4/N3rw0q6T+4Dc6NomNnrunVwt9gmx0bVPavTYK7lF6iNSo2sjowOFRk/9MN4o9Bmh0bWT0YEyo4d+Fm+X+ZDM6NoqMnrmfdwl8imR0bWX0YESo0d+Eu+V+JjE6NosMHribTwg8DmB0bVbXvTAy3hM3oPyoms7owMZHSguet4v4MPinhQXXfsZHSgtetyn8Ahpj0qLricwOpDRgcKip/36HSTsWWHR9QxGBzI6UFb0sN++42Q9LCu6nsLoQEYHyor+8+wBdsl6WFZ0PYXRgYwOFBY967fvMGHPCouuZzA6kNGB0qKH/fodI+1RadH1BEYHioue9ik8QNyT4qJrP6MDGR0oL3rcL+Cj8h6UF13bBUbPu4yHBD4nMLp2S4weeBv3S3xMYnRtZnSgyOiJn8Q7RT4lMrr2yoweeR/3yHxIZnRtZXSg0OiZn8WbhT4jNLp2MjpQavTQD+NtUh+RGl0bxUZPvZIbxD4hNrr2yY0eeyfXyn1AbvTgpV0lePzg6Mlr+17y8MnRoxf3jejRo6Nrj+zo0ffylezBs6Nri/Do2RfzqfCxw6Onr++y9KHTo8cv8IL4keOj63j50ePv5l/5A+dHL1jiXwrGLYiuozVEL7iddw3DNkTXwSqiN1zPLxWjVkTvWOVaLYN2RC9ZZsuYJdE71lkx5OqJ3rDQghFf1UTXcXqix99R/IC/9UTXYYqih19S+Hh/Koqevdbo4f7RFD15scGjfVQVXcfoih57T7GDXdQVPXW5oWN9pix65nojh/pCW/TEBQeO9LWXuonXejl7gL8ULrDu0lfYmqOGuVLjpa+cY+/cXuOlr5hlh4xxq9JLX+v8a+/dXOmlr3X60nubN1/6OvPYu9dWPf1Z2cuXVvx5X2udtP7y5vXRzwjQ3rz+877Weu43fsS+JjziedmHbGvGM56Tfcyupjxkf/Y5m6r/Q+7d5ihzmk+69LXz2GetadRrdmUftqRBn/e11qY8w5qPu/R1/LHP29C0S1+HR5rXfGL0YzMNbD7x8/7qiI/82N1Mfdhaj4UfvZfJj7s/+/CtzH7eWreHB2xk/hPXLd0Z60C88trslGVA3rm+7w7aBOepa63PytOWAHvvpey8FeBevNZ6Lw99PfPZbBP/37u+YXQgowMZHcjoQEYHMjqQ0YGMDmR0IKMDGR3I6EBGBzI6kNGBjA5kdCCjAxkdyOhARgcyOpDRgYwOZHQgowMZHcjoQEYHMjqQ0YGMDmR0IKMDGR3I6EBGBzI6kNGBjA5kdCCjAxkdyOhARgcyOpDRgYwOZHQgowMZHcjoQEYHMjqQ0YGMDmR0IKMDGR3I6EBGBzI6kNGBjA5kdCCjAxkdyOhARgcyOpDRgYwOZHQgowMZHcjoQEYHMjqQ0YGMDmR0IKMDGR3I6EBGBzI6kNGBjA5kdCCjAxkdyOhARgcyOpDRgYwOZHQgowMZHcjoQEYHMjqQ0YGMDmR0IKMDGR3I6EBGBzI6kNGBjA5kdCCjAxkdyOhARgcyOpDRgYwOZHQgowMZHcjoQEYHMjqQ0YGMDmR0IKMDGR3I6EBGBzI6kNGBjA5kdCCjAxkdyOhARgcyOpDRgYwOZHQgowMZHcjoQEYH+h/FOriPguhV6QAAAABJRU5ErkJggg==",
      "text/plain": [
       "<PIL.Image.Image image mode=L size=500x500>"
      ]
     },
     "execution_count": 94,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "t0 = time.time()\n",
    "iterates = np.copy(zz)\n",
    "active = np.ones_like(iterates,dtype=\"bool\")  # Grid of True\n",
    "for _ in range(1000):\n",
    "    iterates[active] = q(iterates[active])  # Apply f to each pixel where we are supposed to proceed\n",
    "    active[active] &= np.abs(iterates[active]) < 2  # In future iterations, only continue on the pixels where abs < 2\n",
    "                                                     # (but only check this for points that remain active)\n",
    "t1 = time.time()\n",
    "print(t1-t0,\"seconds\")\n",
    "PIL.Image.fromarray(255*(1 - active.astype(\"uint8\")))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cbe6745d",
   "metadata": {},
   "source": [
    "### $g(z) = z^2 - 0.52 + 0.53i$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 91,
   "id": "9c6a8ac8",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.4580848217010498 seconds\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAfQAAAH0CAAAAADuvYBWAAAWB0lEQVR4nO2d25rjLK+EpfX893/LWgeAEDvHnpkvhlTVQXfHu6TzWkIIgdWEQtP/vf0BqO+L0AFF6IAidEAROqAIHVCEDihCBxShA4rQAUXogCJ0QBE6oAgdUIQOKEIHFKEDitABReiAInRAETqgCB1QhA4oQgcUoQOK0AFF6IAidEAROqAIHVCEDihCBxShA4rQAUXogCJ0QBE6oAgdUIQOKEIHFKEDitABReiAInRAETqgCB1QhA4oQgcUoQOK0AFF6IAidEAROqAIHVCEDihCBxShA4rQAUXogCJ0QBE6oAgdUIQOKEIHFKEDitABReiAInRAETqgCB1QhA4oQgcUoQOK0AFF6IAidEAROqAIHVCEDihCBxShA4rQAUXogCJ0QBE6oAgdUIQOKEIHFKEDitABReiAInRAETqgCB1QhA4oQgcUoQOK0AFF6IAidEAROqAIHVCEDihCBxShA4rQAUXogCJ0QBE6oAgdUIQOKEIHFKEDitABReiAInRAETqgCB1QhA4oQgcUoQOK0AFF6IAidEAROqAIHVCEDihCBxShA4rQAUXogCJ0QBE6oAgdUIQOKEIHFKEDitABReiAInRAETqgCB1QhA4oQgcUoQOK0AFF6IAidEAROqAIHVCEDihCBxShAwoLur79AfYQFnQjdRE06Gpvf4IthAWdli4iaNDZqosIFnQVMVNiR4KuGbvA2ztUaJNYG9T/PNP/3v4A35SJiqGbuYBBFxVTEXRDh4KuBtacrQQUyCUxeseCnlMzCt+qQ0BHh9wLAnpRbs/h2/Wfh645ZBcVEbOamwE2/5+HboWx5bbcksUDM4fospmYJpP3TjoycgBLFxGN3TTNto+cmQOAbuVXbtHTJjXccA4AepLllrzUUSDXUyC06UXJqRvzsQiWbiZmIipmKibtiIsCRnUI0EsfzUREVdRE812QmaNRPx76PUtVT8h4W64qqgpZUnF8m36n66VWmFv+u4A2UUDq21n6U1er6in11aklXk83SD4hZ+c8UYek7aA/HfG25KilJlz9OvVi4YeZWeceaOkbyCHcpq+ior2jNq1eoF62H2zBc+4bQrdq6veYW5Noay/lw2v5csmv19gv1UnSvb+uxEn9749Hq6gmlz0crsGQg/9ofTteFn6/6N1EMss7NEw0BWOmQ9vcbpmE6ZY2wzXq+0EvHSrTO2Mi5bZI/Np9k5smbrL8Ey8Lvx/0Z1NQUqerRO/XV80npNR7qIBHM/Ttgld15re9rpapKzKc1tlw3GPabQHSbpaeiJfg+hYTK2auVq+Qz196bgPlLSL7WbpIyZY+sfS2rS5njswNMek6ajdLl+ygb1tiaQuKTC+wwoVsc23XT/fxzlk6drKlvztCyD+7cZTzmvaDnjInqnlcpE+X9cB0tq1umBo9/ft+Tdw64vZPG0K1eS+7tNwrowZv2Tds04N6OClQDyjn7PyIRfQOznwj6HMQoeAhday1z65eX2qBHZv6Hm16DK6a3GsdQcvzkXQYIJnh02tTx8u2t9oDus1DsqF9VzPt2nHLd0mHMc9qGYL7/gBI7eHnfJqR51Kzuu5YSKBLm3izMl1J+lMXaTnAwTXXHtAbXG3g3eZUa5JdrfvseUM/irYM4U1h83P7BHJTjfm1JaUuM+ct97RLl7djMt8Fegsrx2widZ0Y01z6ZK0TGKrg+vLIybB6aRxAkW/0j5dhkzpmJiIzH12baEtTFYYJC17yXF6P8T6uaxd5C7qarOOohvqHhlm8odcQuU0CgvakbsoDml7vso3FqLGjXrz7snuVM/WiqeJZpHHnofQ9voGIqBjuBPWX2nQts4Ylt8LNgHgqdHRHL7dmk8dLlguVc8daOb1RYPWzeqlhi+5VQ6o8bdYaj10PnAzyCjh/KbMLmCB799egi8RwLalkZupizbl46kllepvemVPP/h02PfNSmz72oM0s5UZjBs0rXT9dqhZMaMRdn+MwhPdmuKb+WvReU68i4neBdXs+eHdzl5FTsfVibeNgg6Uj6xXoOYjqW9/q1LtG98LW+3EWr2zON432bcjkJDi94d69u6TS22YeOW3SavP+Wppx3OOz7i/ttpYr4o6wibzTZbOagBm7YoNtLiqiRHvmsxzbpMdnHE9/wb0Xx9vkyb2z3vTfwvG94nM5PPZrGvM4VNtQv8wHIuiV5Ezujbe+WONOKbH1gnikFnp0mjI74cSxTTeFXi1S5PttupYKGY2ktUwkT7eDXdeui6mbbojM62pxIjmZm3YN/fxr5L+/SsGX3XvKrlrfTldvrK1ft2ndS3DtZX8I9MsYTLs1vNcqGadW+xM/7Qq+795VykS1ppNefHpaods0DakvrC4g0fYvazvlJaPXnT27qPqBPz/o+lX3riKauloyDHi0ba+a3piAVPNt5r9tnOHSNBHmPzrVRad+nfl3/8FartLky6w2uiafG9Tq3KNvzz7dNPYLLL9r22Nbfbq8/1lkf+Id8lVLD+G61Qypeaje7FpK62pS9d7xiK1YcuzRVeafbypVm+cFFoffP3QffT0jl2Grw1IN25ve9VKlU9a6h3onWHAlbRFdXwefxnM0H21illeVvAXzzmfdUS8VUZj/blMwd7x84dZEBLGXX24CNf/p6n13sutywUmh1eV/oWdm914rlyosHH/TvbpCXi9RfxVwplpDuS6/nza1iMzKqpIpgZDSNjcGevxyJ+Z5dolDWjIiXeo07CyQYuq1DM5JjAiHy07/19J9LInAfKWLCgut/imnHT78b9tpk7r3WRZmbmc1Nd9trJ6jJb16R79FrDTtYavp2h7qBAxdZAx31+vVsEW9Ha7y4x73jz7YSh6ujNxqs7e9rbSsGjgkDG4tWliucuQg7TbQxbpJygl73wTXPzVsCr22evAwuzXeVmVp0XRgDd/y236Ky89eZXKXNj33sZuMue/yv+KUhrpRaiDnEVwHZWzPdSTntVqf0u8eN6yy+JtrH0vXvMZE9y1q73zzwV3w70F82jvi7LaEo+u27DRyG3HjE59IXLYJ5ErmtMTOdXOcElELpGqZazXrlHi5mVWZ+ud802kM1uYfNrw4EPwm0D1ctu5rNDfAYUK6iNSd4YTeyptf6Rxb1N1Vv37ZoEvbZTvOxW/i3k3SsFrzfed8STjsU95cxIso6pWHoMz6Y569Q76qF34cx3wj7xRmKbebyqrQfdJcRMrOUvlcj/RZzLN6ySuy974PLRnkIy19F/fuX2IkVFNzTYC3zIGFzH06Idwz42GLTyFLS/B7ssQQ2gQg52gfS1/AnM0kj3FYm2/th1dG3vMCrKiV7ar3DtP7uls5jvs+lr767j48myUXWNY8SxcGdifeCO59AKdtbVQC3pq69YGAg8DvZOnTDzNr0Md0TRO9L/+l4uxvkF91FcLQq69zdf2mG2r/Dzu9E+ILG7ZNnXNuje925WeX6su64ij+9t9j0FEf1jUFXJv5sSkvnvlRyrx/vkQ//NsM5jy58NvapJ/+ULFNL199RaBhRxgxfbq+/8C8B9v0Lk8agTkP+nICSptwKfOeU631bN2ZW+9V/5ycX1alNR/VPUNn+SUpoyHRv+vcAUuN2kon7CGXcehn3FsbjoO+yJM+q2tRyT5Q8eGbP4M+ZAe6q5/41YnIge59rFXxuU9jebOISB5dWc2Qun6vUEjV3wIfkzz76rjbNSVEljUSs5b3b+GMfcJ2x3E6zdJTJmQ+eioyraz7O+Ze4zy5zKmTmo+z9JpyDSoj3I9WnLsnz71OU3nnfXsiR0Jf2G4d/vinmg7o1p1nfn+Hfujx64/7/+3b+YDafOd5Oq1NF5FJICcibSLun6o+OmZWcHNgs36kpS+mPK067H+rMMtycrsdWDpzJnS56ED9dQdt1usrbzqUyo8VeAfoSPcuIuM37a/tz+aXxfP7y3v6p7MRBnJf1oeFwv4k5ZrrZWyeAJpc+VDoJ1l6GistidEPnbM/hZEG5aSvjK8hXDMjTkQPnMN40J0aq+Em1TIhnppa6gfVy5byp37ZM50EDdMK6+11kKXnZeXc3of9Fo/9k+vnJ7uYhERc0cqeT5ygfhD0tkpiGGtr6lafZubMoquw8G71slaPXX6KM3SUb4rVqCrB+w7PXX1w0VrXOl+epDlyuPiyfndnbVT3fldWf6/gPsjDm6TlKW3RJgSn4a5gOVg7W88q7NlFW32YD2oXg4zLxWnXBj+J5OzjakFdmW3bZZtMxZiXYG/U9p9k6ZZtzCeLxtC6+aq9OdZ7NRSzmulwNZ8/Y937Dih1PTCwD/OjAjkRkbiiYzLz9GivlXOWz0u92Rghdr2zrvIy3BGT4b6xLmv56d7SSZYuItl2rVqdhF8TrcdE4wHx+uXqdYE5kbW/iJG85X67dIZtqqU52kMntelJuUPWDG59+C90MVXVZmd28V8+pj77S5qMjX+oEvH5gvNjd2Kbr/o09y7eulqYYPDh61zt1n4a6yQAzMuIPniHlNxpUvVq9scTLv4DnWfpfXS2/ge0Jlzm33bnclcOoblEtXoZnHu4Z2pGuD6iZJum/UToEg306vPHOOxGrdMSutSOYnz3sExt83n6hRG0vxPe1nHuPenOSp5lFfcmHndNfPbIvD5OoDm6JGtKpjaeNZqzJnvfaLnoE6Hn8NhuzBusBRVu9RXUsKawNI2G5cY5LxtabTmkBIe/UneyXZVWcjS/DfUTobsu8y7ZG6uatIu/pjNna8U2V+5zPd5I+MO6Q46m/SBa3lbiEZtEcXJgP10uk+69vDktr0tWz/rVaaJrd+SlGfdHAeS9fi1Zf5TxGRLbxE/nWXqymepiL0Iv05wYSRtSc6DJQV/l6bo+Wx1vabtm1x3GMQGwja1vc/c90udPXZeY8cEY9e7WJNM+GS/NIbpX0kgzmhszQ8v7rhmgmaWB3tGB0B/0fsbsmE6Gxrp+f1nPqgQFkmui2j7/EMtPZN0B5bXFXy/oPPcuObt1QyXiqqGcio1J8M7RlyVHpTLvLhj+vvTZpWVoinBUvejyrZLKAy390TClfs7V9jmZauiWM/2+zmSfhE+f5lJtU1E/Vfgvvs7gQEt/1OG98USWpjMQnxRWdmh4sGNpXdxWPxGLFdIqomGgtdv1PZ1o6Y9UwrCrf3RM5IcgruRemwZ99liRTlaTxSGH602ONhHid7tzR1r6A9Vil2kD2uXhShMdyyTiucNjnJdKgUebs+sjgpTqS87iq7b329A9D6ar3PdwK1Q/ng8QETWL3fraNC8eECFWHhVu5REjdbSmOfnrvNNb/7h7/5gH64bUrbhbH1HzFiL46tLnmi4rPbxjyRK0Y3HeVfx+Tc2vQ/+kvo+eUPvDXNR/1kNC4L7KB8YGv3b7qqVbOfmdxeJPzL3/M6lYP5wavHvtUdX4vlqlLlPAnoLx1IBvb6OHB2MI/1bglt4vBd7E5Bb+GJuJxRMDvLLCrxcr+ryJmDyx5nuCtvRmmDxtiN2qqH41ymWoEN11jgemxZferr/AHdzSR8Btm15SMYOVL9Und3yFg27spdbUfV+EHl/UIbjWv0vOtfhdcJGTCVct4Zu2Pr/Z84/+jUeChT4rgahWmgstGlIl5rsa7PFHwfmL+JDG2GPrMvhfFSp0HR/5IzHt0ry+H2PndK8Gpz5fwOJV6r+dkVvLvEMV1SZPcr7u4iGdg/o6K1HJVZK1OU/JPb2Ry/2vhAo9pVeHL147G/em9x4fC79TFV4ir55rN48HbDLV8UtCde8iMvbSxzA7evhbhRsJsvfBJz0y7+1F3/9d4Vq6zh/a1S0dWO0yxQBmctEST2Y/TEZ0XrczWOhWyqK67cW9Z36hSKKMmKUh08mg2+RK6bZSlcl6c/bKCJsIsHv3VFmPos/K2eCkY7QWUuqmXlM51EyEi20gUOjzHlvWhNb8oY5hvC0eUmc4ln3Nld8XZu69PPp8qjLG1gdzd7Lt6fDdFw4FtXRZjpItnvxU8ujB5EPjMNRBDB1+qanX94UMfZFCL23zsG8yGOqpt87/Dy5dc072bz/1P9EmH+MtXXr4xe52tmupi1+ZukTPsMliQ9jQl4UQfRK+299cQrrOX93SHr7PV40ZyBUtYjmtf1w93MOPaA24z/LNTntX2NCXE06LVc73t9UxJr0V7zMpeS7YjJyI1CGvcbvvHy3UxhdDzVVz2e16cNjQReauuJ3h0mqWfLfuVU7Xlkz9Po49aZ/o4j3NnPwkvzbuunXx7ZDT0rOWYJoFg4aHPVypPHx9P+a09KKlNZdUTR0Ev3U5X2Bww8d4gUbv88rmsD806mEwbQzUVwotv25HHdS953ommabO+mNr/lzVluM0nVTS06XsSYndlwRn6ertstZCOZGP1VBaz3nSTt+8Rb4rOOgyGU4ZMmrT8dJ4w3x+E9OwZMlmQoMe56k0itNXJ9ONtHS/b0ZmGlqQ7bijtelxodj0ehxb0Uk6Ja9CYtMc3ext6sW3Yw4HPReel9S6587HmvX2tPxSbz9Pt7zDavjlTcFBb/tcodmNZcup4rk5Jxn5bau9Hpt9WXDQW9+uohojrvRjjMCsNPSPLXe7/poAQreyhFspSk/l58N8puE8760/xrgddsjoPRZJSHi8z611njZMqz4VGnQpKwSEmSvO/I5N3greNxeae/cul0l65FpZVnA99+H3dLyr+gtpdui+lvueqZR/Lzz37irTXD48SPsHhebeg3zdz+LZMXy7ILv3HLWPS33+/jeCa+m1IY8Zuh1zKf9cuNClJlwsp1g3evDxfypk6LmPlmN282zsz3PHbdOz+tWkEL4PZEsXqZB/3rqj0KH3QjB0eOjx0dowQocue454/7cCTsNmaRhnAxE89LcenvOm4KGn1hwLO0S/dC3vsUGNtIEHciZSnniLwxzd0pOQrFwE3tKTsJATehJUaobQk8AaOUIXoaVjipZO/bgIHVCEDihCBxShA4rQAUXogCJ0QBE6oAgdUIQOKEIHFKEDitABReiAInRAETqgCB1QhA4oQgcUoQOK0AFF6IAidEAROqAIHVCEDihCBxShA4rQAUXogCJ0QBE6oAgdUIQOKEIHFKEDitABReiAInRAETqgCB1QhA4oQgcUoQOK0AFF6IAidEAROqAIHVCEDihCBxShA4rQAUXogCJ0QBE6oAgdUIQOKEIHFKEDitABReiAInRAETqgCB1QhA4oQgcUoQOK0AFF6IAidEAROqAIHVCEDihCBxShA4rQAUXogCJ0QBE6oAgdUIQOKEIHFKEDitABReiAInRAETqgCB1QhA4oQgcUoQOK0AFF6IAidEAROqAIHVCEDihCBxShA4rQAUXogCJ0QBE6oAgdUIQOKEIHFKEDitABReiAInRAETqgCB1QhA4oQgcUoQOK0AFF6IAidEAROqAIHVCEDihCBxShA4rQAUXogCJ0QBE6oAgdUIQOKEIHFKEDitABReiAInRAETqgCB1QhA4oQgcUoQOK0AFF6IAidEAROqAIHVCEDihCBxShA4rQAUXogCJ0QBE6oP4fLZGFBdwbUJkAAAAASUVORK5CYII=",
      "text/plain": [
       "<PIL.Image.Image image mode=L size=500x500>"
      ]
     },
     "execution_count": 91,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "t0 = time.time()\n",
    "iterates = np.copy(zz)\n",
    "active = np.ones_like(iterates,dtype=\"bool\")  # Grid of True\n",
    "for _ in range(1000):\n",
    "    iterates[active] = g(iterates[active])  # Apply f to each pixel where we are supposed to proceed\n",
    "    active[active] &= np.abs(iterates[active]) < 2  # In future iterations, only continue on the pixels where abs < 2\n",
    "                                                     # (but only check this for points that remain active)\n",
    "t1 = time.time()\n",
    "print(t1-t0,\"seconds\")\n",
    "PIL.Image.fromarray(255*(1 - active.astype(\"uint8\")))"
   ]
  }
 ],
 "metadata": {
  "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.10.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}