{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import pandas as pd\n",
    "import warnings\n",
    "warnings.simplefilter(action='ignore', category=FutureWarning)\n",
    "\n",
    "%config InlineBackend.figure_format ='retina'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Generating Fractal From Random Points - The Chaos Game"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This notebook details how to generate fractals from random points.\n",
    "\n",
    "To get started, checkout this video by [Numberphile](https://www.youtube.com/watch?v=kbKtFN71Lfs)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Initial Definitions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The goal of this exercise is to break down the individual steps of the procedure into functions.\n",
    "\n",
    "To start we generate a list of points that represent our starting point and our verticies."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "def placeStartpoint(npts,fixedpts):\n",
    "    #Start Point\n",
    "    #start = (0.5,0.5)\n",
    "    start = (np.random.random(),np.random.random())\n",
    "\n",
    "    if fixedpts == []: #generates a set of random verticies\n",
    "        for i in range(npts):\n",
    "            randx = np.random.random()\n",
    "            randy = np.random.random()\n",
    "            point = (randx,randy)\n",
    "            fixedpts.append(point)\n",
    "            \n",
    "    return (start,fixedpts)\n",
    "        "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next we build our point chooser which will select a random vertex."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def choosePts(npts,fixedpts,frac):\n",
    "    #chooses a vertex at random\n",
    "    #further rules could be applied here\n",
    "    roll = np.floor(npts*np.random.random())\n",
    "    point = fixedpts[int(roll)]\n",
    "    \n",
    "    return point"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then we iterate on our points."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "def placeItteratePts(npts,itt,start,fixedpts,frac):\n",
    "    ittpts = []\n",
    "    \n",
    "    for i in range(itt):\n",
    "        point = choosePts(npts,fixedpts,frac) #chooses a vertex at random\n",
    "#        halfway = ((point[0]+start[0])*frac,(point[1]+start[1])*frac) #calculates the halfway point between the starting point and the vertex\n",
    "        halfway = ((point[0]-start[0])*(1.0 - frac)+start[0],(point[1]-start[1])*(1.0 - frac)+start[1])\n",
    "        ittpts.append(halfway)\n",
    "        start = halfway #sets the starting point to the new point\n",
    "        \n",
    "    return ittpts"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Lastly, we need to plot all of our points."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "def plotFractal(start,fixedpts,ittpts):\n",
    "    # set axes range\n",
    "    plt.xlim(-0.05,1.05)\n",
    "    plt.ylim(-0.05,1.05)\n",
    "    plt.axis('equal')\n",
    "    \n",
    "    #plots the verticies\n",
    "    plt.scatter(np.transpose(fixedpts)[0],np.transpose(fixedpts)[1],alpha=0.8, c='black', edgecolors='none', s=30)\n",
    "    #plots the starting point\n",
    "    plt.scatter(start[0],start[1],alpha=0.8, c='red', edgecolors='none', s=30)    \n",
    "    #plots the itterated points\n",
    "    plt.scatter(np.transpose(ittpts)[0],np.transpose(ittpts)[1],alpha=0.5, c='blue', edgecolors='none', s=2)\n",
    "    \n",
    "    plt.show()\n",
    "    return"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Generate Fractal"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Once we have all of our functions, we can build a central function that calls everything we might need to generate an image."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "def GenerateFractal(npts,frac,itt,reg=False):\n",
    "    #Error Control\n",
    "    if npts < 1 or frac >= 1.0 or frac <= 0.0 or type(npts) is not int or type(frac) is not float or type(itt) is not int:\n",
    "        print(\"number of points must be a positive integer, compression fraction must be a positive float less than 1.0, itt must be a positive integer\")\n",
    "        return\n",
    "    if frac > 0.5:\n",
    "        print(\"Warning: compression fractions over 1/2 do not lead to fractals\")\n",
    "\n",
    "    #Initilize Verticies\n",
    "    if not reg:\n",
    "        fixedpts = [] #Random Verticies\n",
    "    else:\n",
    "        if npts == 3:\n",
    "            fixedpts = [(0.0,0.0),(1.0,0.0),(0.5,0.5*np.sqrt(3.0))] #Equilateral Triangle (npts = 3)\n",
    "        elif npts == 4:\n",
    "            fixedpts = [(0.0,0.0),(1.0,0.0),(1.0,1.0),(0.0,1.0)] #Square\n",
    "        elif npts == 5:\n",
    "            fixedpts = [(0.0,2./(1+np.sqrt(5.))),(0.5-2./(5+np.sqrt(5.)),0.0),(0.5,1.0),(0.5+2./(5+np.sqrt(5.)),0.0),(1.0,2./(1+np.sqrt(5.)))] #Regular Pentagon\n",
    "        elif npts == 6:\n",
    "            fixedpts = [(0.0,0.5),(1./4,0.5+.25*np.sqrt(3.)),(3./4,0.5+.25*np.sqrt(3.)),(1.0,0.5),(3./4,0.5-.25*np.sqrt(3.)),(1./4,0.5-.25*np.sqrt(3.))] #Regular Hexagon\n",
    "        elif npts == 2:\n",
    "            fixedpts = [(0.0,0.0),(1.0,1.0)] #Line\n",
    "        elif npts == 1:\n",
    "            fixedpts = [(0.5,0.5)] #Line\n",
    "        elif npts == 8:\n",
    "            fixedpts = [(0.0,0.0),(1.0,0.0),(1.0,1.0),(0.0,1.0),(0.5,0.0),(1.0,0.5),(0.5,1.0),(0.0,0.5)] #Carpet\n",
    "        else:\n",
    "            print(\"No regular polygon stored with that many verticies, switching to default with randomly assigned verticies\")\n",
    "            fixedpts = [] #Random Verticies\n",
    "\n",
    "    #Compression Fraction\n",
    "#    frac = 1.0/2.0 #Sierpinski's Triangle (npts = 3)\n",
    "#    frac = 1.0/2.0 #Sierpinski's \"Square\" (filled square, npts = 4)\n",
    "#    frac = 1.0/3.0 #Sierpinski's Pentagon (npts = 5)\n",
    "#    frac = 3.0/8.0 #Sierpinski's Hexagon (npts = 6)\n",
    "\n",
    "        \n",
    "    if len(fixedpts) != npts and len(fixedpts) != 0:\n",
    "        print(\"The number of verticies don't match the length of the list of verticies. If you want the verticies generated at random, set fixedpts to []\")\n",
    "        return\n",
    "    if len(fixedpts) != 0:\n",
    "        print(\"Fractal Dimension = {}\".format(-np.log(npts)/np.log(frac)))\n",
    "    \n",
    "        \n",
    "    (start, fixedpts) = placeStartpoint(npts,fixedpts)\n",
    "    \n",
    "    ittpts = placeItteratePts(npts,itt,start,fixedpts,frac)\n",
    "    \n",
    "    plotFractal(start,fixedpts,ittpts)\n",
    "    \n",
    "    return\n",
    "    \n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Make A Fractal"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "# Call the GenerateFractal function with a number of verticies, a number of itterations, and the compression fraction\n",
    "# The starting verticies are random by default. An optional input of True will set the verticies to those of a regular polygon.\n",
    "GenerateFractal(3,.5,5000)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Regular Polygons"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There are some famous fractals that can be generated this way that are rooted in regular polygons. These are often named after the Polish mathematician Waclaw Sierpinski but have existed long before he categorized them. For example [Sierpinski's Triangle](https://en.wikipedia.org/wiki/Sierpi%C5%84ski_triangle) and [Sierpinski's Carpet](https://en.wikipedia.org/wiki/Sierpi%C5%84ski_carpet)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "GenerateFractal(3,.5,50000,True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "GenerateFractal(5,1./3,50000,True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    },
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "GenerateFractal(6,3./8,50000,True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "GenerateFractal(8,1./3.,50000,True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true
   },
   "source": [
    "#### Exploring Further: Dimension"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Fractals have, by definition, [fractional dimension](https://en.wikipedia.org/wiki/Fractal_dimension). We can model what that means here.\n",
    "\n",
    "Roughly, fractional dimension is a measure of how much space a fractal fills and in our case will fall between 1 (a line) and 2 (a plane)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "hidden": true,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "GenerateFractal(1,.5,50000,True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "hidden": true,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "GenerateFractal(2,.5,50000,True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "hidden": true,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "GenerateFractal(4,.5,50000,True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true
   },
   "source": [
    "#### Randomness on Large Scales"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To explore the nature of randomness, we can use our same function but increase the number of verticies that are randomly generated to a significantly large number.\n",
    "\n",
    "This is an example of a [Monte Carlo method](https://en.wikipedia.org/wiki/Monte_Carlo_method)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "hidden": true,
    "jupyter": {
     "outputs_hidden": false
    },
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "GenerateFractal(10,.5,100)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "hidden": true,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "GenerateFractal(10,.5,5000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "hidden": true,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "GenerateFractal(100,.5,5000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "hidden": true,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "GenerateFractal(100,.5,100000)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Learn More:\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[Chaos Game Wiki](https://en.wikipedia.org/wiki/Chaos_game)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[Numberphile Video](https://www.youtube.com/watch?v=kbKtFN71Lfs)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[Chaos in the Classroom](http://math.bu.edu/DYSYS/chaos-game/chaos-game.html)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[Chaos Rules!](http://www.maa.org/sites/default/files/pdf/upload_library/2/Devaney%202005.pdf)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[Barnsley Fern](https://en.wikipedia.org/wiki/Barnsley_fern)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Modeling Life"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There are many examples of fractals in nature. Here we have generated a fern. If we use the specific values below this is known as [Barnsley's Fern](https://en.wikipedia.org/wiki/Barnsley_fern)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def makeFern(f,itt):            \n",
    "    colname = [\"percent\",\"a\",\"b\",\"c\",\"d\",\"e\",\"f\"]\n",
    "    print(pd.DataFrame(data=np.array(f), columns = colname))\n",
    "    \n",
    "    x,y = {0.5,0.0}\n",
    "    xypts=[]\n",
    "    if abs(sum(f[j][0] for j in range(len(f)))-1.0) < 10^-10:\n",
    "        print(\"Probabilities must sum to 1\")\n",
    "        return\n",
    "    for i in range(itt):\n",
    "        rand = (np.random.random())\n",
    "        cond = 0.0\n",
    "        for j in range(len(f)):\n",
    "            if  (cond <= rand) and (rand <= (cond+f[j][0])):\n",
    "                x = f[j][1]*x+f[j][2]*y+f[j][5]\n",
    "                y = f[j][3]*x+f[j][4]*y+f[j][6]\n",
    "                xypts.append((x,y))\n",
    "            cond = cond + f[j][0]\n",
    "            \n",
    "    xmax,ymax = max(abs(np.transpose(xypts)[0])),max(abs(np.transpose(xypts)[1]))\n",
    "    plt.axes().set_aspect('equal')\n",
    "    color = np.transpose([[abs(r)/xmax for r in np.transpose(xypts)[0]],[abs(g)/ymax for g in np.transpose(xypts)[1]],[b/itt for b in range(itt)]])\n",
    "    \n",
    "    plt.scatter(np.transpose(xypts)[0],np.transpose(xypts)[1],alpha=0.5, facecolors=color, edgecolors='none', s=1)\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### For Barnsley's Fern:\n",
    "Use the following values\n",
    "\n",
    "|__Percent__|__A__|__B__|__C__|__D__|__E__|__F__|\n",
    "|---|---|---|---|---|---|---|\n",
    "|0.01|0.0|0.0|0.0|0.16|0.0|0.0|\n",
    "|0.85|0.85|0.04|-0.04|0.85|0.0|1.60|\n",
    "|0.07|0.20|-0.26|0.23|0.22|0.0|1.60|\n",
    "|0.07|-0.15|0.28|0.26|0.24|0.0|0.44|\n",
    "\n",
    "Of course, this is only one solution so try as changing the values. Some values modify the curl, some change the thickness, others completely rearrange the structure."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "f = ((0.01,0.0,0.0,0.0,0.16,0.0,0.0),\n",
    "     (0.85,0.85,0.08,-0.08,0.85,0.0,1.60),\n",
    "     (0.07,0.20,-0.26,0.23,0.22,0.0,1.60),\n",
    "     (0.07,-0.15,0.28,0.26,0.24,0.0,0.44))\n",
    "\n",
    "makeFern(f,5000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "jupytext": {
   "formats": "ipynb,md"
  },
  "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.8.10"
  },
  "nav_menu": {},
  "toc": {
   "navigate_menu": true,
   "number_sections": false,
   "sideBar": false,
   "threshold": 6,
   "toc_cell": true,
   "toc_section_display": "block",
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}