{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Visualizing Attractors\n",
    "\n",
    "An [attractor](https://en.wikipedia.org/wiki/Attractor#Strange_attractor) is a set of values to which a numerical system tends to evolve. An attractor is called a [strange attractor](https://en.wikipedia.org/wiki/Attractor#Strange_attractor) if the resulting pattern has a fractal structure. This notebook shows how to calculate and plot two-dimensional attractors of a variety of types, using code and parameters primarily from [Lázaro Alonso](https://lazarusa.github.io/Webpage/codepython2.html), [François Pacull](https://aetperf.github.io/2018/08/29/Plotting-Hopalong-attractor-with-Datashader-and-Numba.html), [Jason Rampe](https://softologyblog.wordpress.com/2017/03/04/2d-strange-attractors), [Paul Bourke](http://paulbourke.net/fractals/), and [James A. Bednar](http://github.io/jbednar).\n",
    "\n",
    "\n",
    "## Clifford Attractors\n",
    "\n",
    "For example, a [Clifford Attractor](http://paulbourke.net/fractals/clifford) is a strange attractor defined by two iterative equations that determine the _x,y_ locations of discrete steps in the path of a particle across a 2D space, given a starting point _(x0,y0)_ and the values of four parameters _(a,b,c,d)_:\n",
    "\n",
    "\\begin{equation}\n",
    "x_{n +1} = \\sin(a y_{n}) + c \\cos(a x_{n})\\\\\n",
    "y_{n +1} = \\sin(b x_{n}) + d \\cos(b y_{n})\n",
    "\\end{equation}\n",
    "\n",
    "At each time step, the equations define the location for the following time step, and the accumulated locations show the areas of the 2D plane most commonly visited by the imaginary particle.  \n",
    "\n",
    "It's easy to calculate these values in Python using [Numba](http://numba.pydata.org). First, we define the iterative attractor equation:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np, pandas as pd, datashader as ds\n",
    "from datashader import transfer_functions as tf\n",
    "from datashader.colors import inferno, viridis\n",
    "from numba import jit\n",
    "from math import sin, cos, sqrt, fabs\n",
    "\n",
    "@jit(nopython=True)\n",
    "def Clifford(x, y, a, b, c, d, *o):\n",
    "    return sin(a * y) + c * cos(a * x), \\\n",
    "           sin(b * x) + d * cos(b * y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We then evaluate this equation 10 million times, creating a set of _x,y_ coordinates visited. The `@jit` here and above is optional, but it makes the code 50x faster."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "n=10000000\n",
    "\n",
    "@jit(nopython=True)\n",
    "def trajectory_coords(fn, x0, y0, a, b=0, c=0, d=0, e=0, f=0, n=n):\n",
    "    x, y = np.zeros(n), np.zeros(n)\n",
    "    x[0], y[0] = x0, y0\n",
    "    for i in np.arange(n-1):\n",
    "        x[i+1], y[i+1] = fn(x[i], y[i], a, b, c, d, e, f)\n",
    "    return x,y\n",
    "\n",
    "def trajectory(fn, x0, y0, a, b=0, c=0, d=0, e=0, f=0, n=n):\n",
    "    x, y = trajectory_coords(fn, x0, y0, a, b, c, d, e, f, n)\n",
    "    return pd.DataFrame(dict(x=x,y=y))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "df = trajectory(Clifford, 0, 0, -1.3, -1.3, -1.8, -1.9)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df.tail()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can now aggregate these 10,000,000 continuous coordinates into a discrete 2D rectangular grid with [Datashader](http://datashader.org), counting each time a point fell into that grid cell:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "\n",
    "cvs = ds.Canvas(plot_width = 700, plot_height = 700)\n",
    "agg = cvs.points(df, 'x', 'y')\n",
    "print(agg.values[190:195,190:195],\"\\n\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A small portion of that grid is shown above, but it's difficult to see the grid's structure from the numerical values.  To see the entire array at once, we can turn each grid cell into a pixel, using a greyscale value from white to black:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ds.transfer_functions.Image.border=0\n",
    "\n",
    "tf.shade(agg, cmap = [\"white\", \"black\"])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As you can see, the most-visited areas of the plane have an interesting structure for this set of parameters. To explore further, let's wrap up the above aggregation and shading commands into a function so we can apply them more easily:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def dsplot(fn, vals, n=n, cmap=viridis, label=True):\n",
    "    \"\"\"Return a Datashader image by collecting `n` trajectory points for the given attractor `fn`\"\"\"\n",
    "    lab = (\"{}, \"*(len(vals)-1)+\" {}\").format(*vals) if label else None\n",
    "    df  = trajectory(fn, *vals, n=n)\n",
    "    cvs = ds.Canvas(plot_width = 300, plot_height = 300)\n",
    "    agg = cvs.points(df, 'x', 'y')\n",
    "    img = tf.shade(agg, cmap=cmap, name=lab)\n",
    "    return img"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And let's load some colormaps that we can use for subsequent plots:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from colorcet import palette\n",
    "palette[\"viridis\"]=viridis\n",
    "palette[\"inferno\"]=inferno"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can now use these colormaps with a pre-selected set of Clifford attractor parameter values (stored in a separate [YAML-format text file](https://raw.githubusercontent.com/pyviz-topics/examples/master/attractors/strange_attractors.yml)) to show a wide variety of trajectories that these equations can form:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import yaml\n",
    "vals = yaml.load(open(\"strange_attractors.yml\",\"r\"), Loader=yaml.FullLoader)\n",
    "\n",
    "def args(name):\n",
    "    \"\"\"Return a list of available argument lists for the given type of attractor\"\"\"\n",
    "    return [v[1:] for v in vals if v[0]==name]  \n",
    "\n",
    "def plot(fn, vals=None, **kw):\n",
    "    \"\"\"Plot the given attractor `fn` once per provided set of arguments.\"\"\"\n",
    "    vargs=args(fn.__name__) if vals is None else vals\n",
    "    return tf.Images(*[dsplot(fn, v[1:], cmap=palette[v[0]][::-1], **kw) for v in vargs]).cols(4)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot(Clifford)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here the values shown are the arguments for the first call to `Clifford(x, y, a, b, c, d)`, with each subsequent call using the _x,y_ location of the previous call.  \n",
    "\n",
    "Randomly sampling the parameter space typically yields much less dramatic patterns, such as all trajectory locations being on a small number of points:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy.random\n",
    "numpy.random.seed(21)\n",
    "num = 4\n",
    "\n",
    "rvals=np.c_[np.zeros((num,2)), numpy.random.random((num,4))*4-2]\n",
    "plot(Clifford, vals=[[\"kbc\"]+list(rvals[i]) for i in range(len(rvals))], label=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If you wish, Datashader could easily be used to filter out such uninteresting examples, by applying a criterion to the aggregate array before shading and showing only those that remain (e.g. rejecting those where 80% of the pixel bins are empty).\n",
    "\n",
    "\n",
    "## De Jong attractors\n",
    "\n",
    "A variety of other sets of attractor equations have been proposed, such as these from [Peter de Jong](http://paulbourke.net/fractals/peterdejong):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "@jit(nopython=True)\n",
    "def De_Jong(x, y, a, b, c, d, *o):\n",
    "    return sin(a * y) - cos(b * x), \\\n",
    "           sin(c * x) - cos(d * y)\n",
    "\n",
    "plot(De_Jong)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Svensson attractors\n",
    "\n",
    "From [Johnny Svensson](http://paulbourke.net/fractals/peterdejong/):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "@jit(nopython=True)\n",
    "def Svensson(x, y, a, b, c, d, *o):\n",
    "    return d * sin(a * x) - sin(b * y), \\\n",
    "           c * cos(a * x) + cos(b * y)\n",
    "\n",
    "plot(Svensson)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Bedhead Attractor\n",
    "\n",
    "From [Ivan Emrich](https://www.deviantart.com/jaguarfacedman) and [Jason Rampe](https://softologyblog.wordpress.com/2017/03/04/2d-strange-attractors):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "@jit(nopython=True)\n",
    "def Bedhead(x, y, a, b, *o):\n",
    "    return sin(x*y/b)*y + cos(a*x-y), \\\n",
    "           x + sin(y)/b\n",
    "\n",
    "plot(Bedhead)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Fractal Dream Attractor\n",
    "\n",
    "From Clifford A. Pickover's book “Chaos In Wonderland”, with parameters from [Jason Rampe](https://softologyblog.wordpress.com/2017/03/04/2d-strange-attractors):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "@jit(nopython=True)\n",
    "def Fractal_Dream(x, y, a, b, c, d, *o):\n",
    "    return sin(y*b)+c*sin(x*b), \\\n",
    "           sin(x*a)+d*sin(y*a)\n",
    "\n",
    "plot(Fractal_Dream)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Hopalong attractors\n",
    "\n",
    "From Barry Martin, here with code for two variants from [François Pacull](https://aetperf.github.io/2018/08/29/Plotting-Hopalong-attractor-with-Datashader-and-Numba.html):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "@jit(nopython=True)\n",
    "def Hopalong1(x, y, a, b, c, *o):\n",
    "    return y - sqrt(fabs(b * x - c)) * np.sign(x), \\\n",
    "           a - x\n",
    "@jit(nopython=True)\n",
    "def Hopalong2(x, y, a, b, c, *o):\n",
    "    return y - 1.0 - sqrt(fabs(b * x - 1.0 - c)) * np.sign(x - 1.0), \\\n",
    "           a - x - 1.0\n",
    "\n",
    "plot(Hopalong1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot(Hopalong2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##  Gumowski-Mira Attractor\n",
    "\n",
    "From [I. Gumowski and C. Mira](http://kgdawiec.bplaced.net/badania/pdf/cacs_2010.pdf), with code and parameters from [Jason Rampe](https://softologyblog.wordpress.com/2017/03/04/2d-strange-attractors) and [Lázaro Alonso](https://lazarusa.github.io/Webpage/codepython2.html):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "@jit(nopython=True)\n",
    "def G(x, mu):\n",
    "    return mu * x + 2 * (1 - mu) * x**2 / (1.0 + x**2)\n",
    "\n",
    "@jit(nopython=True)\n",
    "def Gumowski_Mira(x, y, a, b, mu, *o):\n",
    "    xn = y + a*(1 - b*y**2)*y  +  G(x, mu)\n",
    "    yn = -x + G(xn, mu)\n",
    "    return xn, yn\n",
    "\n",
    "plot(Gumowski_Mira)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##  Symmetric Icon Attractor\n",
    "\n",
    "The Hopalong and Gumowski-Mira equations often result in symmetric patterns, but a different approach is to *force* the patterns to be symmetric, which is often pleasing. Examples from “Symmetry in Chaos” by Michael Field and Martin Golubitsky, with code and parameters from [Jason Rampe](https://softologyblog.wordpress.com/2017/03/04/2d-strange-attractors):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "@jit(nopython=True)\n",
    "def Symmetric_Icon(x, y, a, b, g, om, l, d, *o):\n",
    "    zzbar = x*x + y*y\n",
    "    p = a*zzbar + l\n",
    "    zreal, zimag = x, y\n",
    "    \n",
    "    for i in range(1, d-1):\n",
    "        za, zb = zreal * x - zimag * y, zimag * x + zreal * y\n",
    "        zreal, zimag = za, zb\n",
    "    \n",
    "    zn = x*zreal - y*zimag\n",
    "    p += b*zn\n",
    "    \n",
    "    return p*x + g*zreal - om*y, \\\n",
    "           p*y - g*zimag + om*x\n",
    "\n",
    "plot(Symmetric_Icon)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Interactive plotting\n",
    "\n",
    "If you are running a live Python process, you can use Datashader with [hvPlot](https://hvplot.holoviz.org) to zoom in and see the individual steps in any of these calculations using a [Bokeh](https://bokeh.org) plot:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import hvplot.pandas # noqa\n",
    "\n",
    "df = trajectory(Clifford, *(args(\"Clifford\")[5][1:]))\n",
    "opts = dict(rasterize=True, dynspread=True, cnorm='eq_hist', width=400, height=400, colorbar=False)\n",
    "\n",
    "df.hvplot.scatter(x='x', y='y', **opts)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Each time you zoom in in a live process, the data will be reaggregated, which should take a small fraction of a second for 10 million points.  Eventually, once you zoom in enough you should see individual data points, as we are not connecting the points into a trajectory here. \n",
    "\n",
    "You can also try \"connecting the dots\", which will reveal how the particle jumps discretely from one region of the space to another:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df.hvplot.line(x='x', y='y', **opts)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Again, if you zoom in on a live server, the plot will update so that you can see the individual traces involved. \n",
    "\n",
    "On the live server, you can also explore to find your own parameter values that generate interesting patterns. Here we've made a tiny [Panel](https://panel.holoviz.org) app to regenerate the values given each parameter value, then plot the result:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def hv_clif(a,b,c,d,x0=0,y0=0,n=n):\n",
    "    df = trajectory(Clifford, x0, y0, a, b, c, d, n)\n",
    "    return df.hvplot.scatter(x='x', y='y', **opts)\n",
    "\n",
    "import panel as pn\n",
    "init = zip(['a', 'b', 'c', 'd'], args(\"Clifford\")[6][-4:])\n",
    "widgets = [pn.widgets.FloatSlider(name=n, start=-2.0, end=2.0, value=v) for n,v in init]\n",
    "pn.Row(pn.bind(hv_clif, *widgets), pn.Column(*widgets))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Although many of the regions of this four-dimensional parameter space generate uninteresting trajectories such as single points, you can find interesting regions by starting with one of the _a,b,c,d_ tuples of values in previous Clifford plots, then click on one slider and use the left and right arrow keys to see how the plot changes as that parameter changes. See also this more ambitious [Panel](http://panel.pyviz.org)-based [attractor dashboard](https://examples.pyviz.org/attractors/attractors_panel.html)."
   ]
  }
 ],
 "metadata": {
  "language_info": {
   "name": "python",
   "pygments_lexer": "ipython3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}