{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Iterated function system  with k3d visualization\n",
    "\n",
    "Starting from an array of [0, 0, 0] points, we run a number of iterations that transform each point with one of four functions, with given probabilities (equal by default)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import  numpy as np\n",
    "\n",
    "N = 30000\n",
    "x = np.zeros(N)\n",
    "y = np.zeros(N)\n",
    "z = np.zeros(N)\n",
    "\n",
    "x1 = np.empty_like(x)\n",
    "y1 = np.empty_like(y)\n",
    "z1 = np.empty_like(z)\n",
    "\n",
    "# Sierpinski triangle iterative functions\n",
    "\n",
    "def f1(x,y,z,x1,y1,z1,c):\n",
    "    x1[c] = 1.0/2.0*x[c]\n",
    "    y1[c] = 1.0/2.0*y[c]\n",
    "    z1[c] = 1.0/2.0*z[c]\n",
    "\n",
    "def f2(x,y,z,x1,y1,z1,c):\n",
    "    x1[c] = 1.0/2.0*x[c] + 1/2.0\n",
    "    y1[c] = 1.0/2.0*y[c]\n",
    "    z1[c] = 1.0/2.0*z[c]\n",
    "\n",
    "\n",
    "def f3(x,y,z,x1,y1,z1,c):\n",
    "    x1[c] = 1.0/2.0*x[c] + 1/4.\n",
    "    y1[c] = 1.0/2.0*y[c] + np.sqrt(3)/4\n",
    "    z1[c] = 1.0/2.0*z[c]\n",
    "\n",
    "def f4(x,y,z,x1,y1,z1,c):\n",
    "    x1[c] = 1.0/2.0*x[c] + 1/4.\n",
    "    y1[c] = 1.0/2.0*y[c] + 1./4\n",
    "    z1[c] = 1.0/2.0*z[c] + np.sqrt(3)/4\n",
    "\n",
    "functions = [f1, f2, f3, f4] \n",
    "probabilities = [1/4.]*4\n",
    "assert(len(functions) == len(probabilities))\n",
    "\n",
    "X,Y,Z = x,y,z\n",
    "for i in range(20):\n",
    "    # pick indices for each function to be applied\n",
    "    r = np.random.choice(len(probabilities), size=N, p=probabilities)\n",
    "    for i, f in enumerate(functions):\n",
    "        f(x, y, z, x1, y1, z1, r==i)\n",
    "    x,x1 = x1,x\n",
    "    y,y1 = y1,y\n",
    "    z,z1 = z1,z\n",
    "    if i > 0:\n",
    "        X, Y, Z = np.hstack([X,x]), np.hstack([Y,y]), np.hstack([Z,z])\n",
    "\n",
    "# how much memory are we using, how many points there are\n",
    "print(3*X.nbytes//1024**2,\"MB\",X.shape[0])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we turn separate coordinate array into triplets.\n",
    "Using `numpy` for what could also be written as `list(zip(X, Y, Z))`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "positions = np.vstack([X, Y, Z]).T.copy()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import k3d\n",
    "\n",
    "plot = k3d.plot()\n",
    "point_plot = k3d.points(positions.astype(np.float32), color=0xff0000, point_size=0.003, shader='3d')\n",
    "plot += point_plot\n",
    "plot.display()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Zoom in a little:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot.camera = [-0.59265772150826,\n",
    " 1.0966590944867525,\n",
    " 0.15381683182413644,\n",
    " 0.35173312413637553,\n",
    " 0.35752558265043016,\n",
    " 0.3151305910837551,\n",
    " -0.5602813338387698,\n",
    " -0.7160643522547137,\n",
    " -0.41633720753942915]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.6.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}