chaos_game.ipynb

chaos_game.ipynb

{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "colab": {
      "provenance": [],
      "include_colab_link": true
    },
    "kernelspec": {
      "name": "python3",
      "display_name": "Python 3"
    },
    "language_info": {
      "name": "python"
    }
  },
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "view-in-github",
        "colab_type": "text"
      },
      "source": [
        "<a href=\"https://colab.research.google.com/gist/tueda/5afb14c54400abb669e9d26bb799a85a/chaos_game.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "e5WFZZfEGLww",
        "cellView": "form"
      },
      "outputs": [],
      "source": [
        "# @title カオスゲーム\n",
        "%matplotlib inline\n",
        "%config InlineBackend.figure_format = 'retina'\n",
        "import matplotlib.colors\n",
        "import matplotlib.pyplot as plt\n",
        "import matplotlib.typing as mplt\n",
        "import numpy as np\n",
        "import numpy.typing as npt\n",
        "import inspect\n",
        "from collections.abc import Sequence\n",
        "from numpy import cos, cosh, e, exp, log, log10, sin, sinh, sqrt, pi, tan, tanh\n",
        "from numpy.linalg import inv\n",
        "from typing import Callable, TypeAlias\n",
        "\n",
        "np.set_printoptions(sign=' ', suppress=True)\n",
        "\n",
        "# 黄金比\n",
        "phi = float((1 + sqrt(5)) / 2)\n",
        "\n",
        "\n",
        "def Vec(x: float, y: float) -> np.ndarray:\n",
        "    \"\"\"ベクトルを返す。\n",
        "\n",
        "    Args:\n",
        "        x (float): xの値。\n",
        "        y (float): yの値。\n",
        "\n",
        "    Returns:\n",
        "        numpy.ndarray: ベクトル。\n",
        "    \"\"\"\n",
        "    return np.array((x, y, 1), dtype=np.float64)\n",
        "\n",
        "\n",
        "def I() -> np.ndarray:\n",
        "    \"\"\"恒等行列（identity matrix）を返す。\n",
        "\n",
        "    Returns:\n",
        "        numpy.ndarray: 恒等行列。\n",
        "    \"\"\"\n",
        "    return np.eye(3, dtype=np.float64)\n",
        "\n",
        "\n",
        "def R(theta: float) -> np.ndarray:\n",
        "    \"\"\"回転行列（rotation matrix）を返す。\n",
        "\n",
        "    Args:\n",
        "        theta (float): 回転角。\n",
        "\n",
        "    Returns:\n",
        "        numpy.ndarray: 回転行列。\n",
        "    \"\"\"\n",
        "    c = np.cos(theta)\n",
        "    s = np.sin(theta)\n",
        "    return np.array(((c, -s, 0), (s, c, 0), (0, 0, 1)), dtype=np.float64)\n",
        "\n",
        "\n",
        "def S(sx: float, sy: float) -> np.ndarray:\n",
        "    \"\"\"拡大縮小行列（scaling matrix）を返す。\n",
        "\n",
        "    Args:\n",
        "        sx (float): x軸方向のスケーリング係数。\n",
        "        sy (float): y軸方向のスケーリング係数。\n",
        "\n",
        "    Returns:\n",
        "        numpy.ndarray: 拡大縮小行列。\n",
        "    \"\"\"\n",
        "    return np.array(((sx, 0, 0), (0, sy, 0), (0, 0, 1)), dtype=np.float64)\n",
        "\n",
        "\n",
        "def T(tx: float, ty: float) -> np.ndarray:\n",
        "    \"\"\"平行移動行列（translation matrix）を返す。\n",
        "\n",
        "    Args:\n",
        "        tx (float): x軸方向への変位。\n",
        "        ty (float): y軸方向への変位。\n",
        "\n",
        "    Returns:\n",
        "        numpy.ndarray: 平行移動行列。\n",
        "    \"\"\"\n",
        "    return np.array(((1, 0, tx), (0, 1, ty), (0, 0, 1)), dtype=np.float64)\n",
        "\n",
        "\n",
        "def Sh(shx: float, shy: float) -> np.ndarray:\n",
        "    \"\"\"剪断行列（shear matrix）を返す。\n",
        "\n",
        "    Args:\n",
        "        shx (float): x軸方向の剪断因子。\n",
        "        shy (float): y軸方向の剪断因子。\n",
        "\n",
        "    Returns:\n",
        "        numpy.ndarray: 平行移動行列。\n",
        "    \"\"\"\n",
        "    return np.array(((1, shx, 0), (shy, 1, 0), (0, 0, 1)), dtype=np.float64)\n",
        "\n",
        "\n",
        "def A(\n",
        "    a11: float, a12: float, a13: float, a21: float, a22: float, a23: float\n",
        ") -> np.ndarray:\n",
        "    \"\"\"アフィン行列（affine matrix）を返す。\n",
        "\n",
        "    Args:\n",
        "        a11 (float): (1, 1)成分。\n",
        "        a12 (float): (1, 2)成分。\n",
        "        a13 (float): (1, 3)成分。\n",
        "        a21 (float): (2, 1)成分。\n",
        "        a22 (float): (2, 2)成分。\n",
        "        a23 (float): (2, 3)成分。\n",
        "\n",
        "    Returns:\n",
        "        numpy.ndarray: アフィン行列。\n",
        "    \"\"\"\n",
        "    return np.array(((a11, a12, a13), (a21, a22, a23), (0, 0, 1)), dtype=np.float64)\n",
        "\n",
        "\n",
        "def _is_real_scalar(x: object) -> bool:\n",
        "    return isinstance(x, (int, float, np.integer, np.floating)) and not isinstance(\n",
        "        x, bool\n",
        "    )\n",
        "\n",
        "\n",
        "def _is_affine_like(x: object) -> bool:\n",
        "    try:\n",
        "        a = np.asarray(x)\n",
        "        return a.shape == (3, 3) and (\n",
        "            np.issubdtype(a.dtype, np.integer) or np.issubdtype(a.dtype, np.floating)\n",
        "        )\n",
        "    except Exception:\n",
        "        return False\n",
        "\n",
        "\n",
        "def _is_sequence_but_not_ndarray(x: object) -> bool:\n",
        "    return isinstance(x, Sequence) and not isinstance(\n",
        "        x, (np.ndarray, str, bytes, bytearray)\n",
        "    )\n",
        "\n",
        "\n",
        "def _parse_maps(\n",
        "    x: object, weights: list[float], maps: list[npt.NDArray[np.float64]]\n",
        ") -> bool:\n",
        "    match x:\n",
        "        case m if _is_affine_like(m):\n",
        "            weights.append(1.0)\n",
        "            maps.append(np.asarray(m, dtype=np.float64))\n",
        "        case (m, w) if _is_affine_like(m) and _is_real_scalar(w):\n",
        "            weights.append(float(w))\n",
        "            maps.append(np.asarray(m, dtype=np.float64))\n",
        "        case (w, m) if _is_affine_like(m) and _is_real_scalar(w):\n",
        "            weights.append(float(w))\n",
        "            maps.append(np.asarray(m, dtype=np.float64))\n",
        "        case _:\n",
        "            return False\n",
        "    return True\n",
        "\n",
        "\n",
        "AffineMapLike: TypeAlias = (\n",
        "    npt.ArrayLike | tuple[npt.ArrayLike, float] | tuple[float, npt.ArrayLike]\n",
        ")\n",
        "\n",
        "# (x: float, y: float) -> float\n",
        "# (sample: int, x: float, y: float) -> float\n",
        "# (iteration: int, sample: int, x: float, y: float) -> float\n",
        "ColorFunc2: TypeAlias = Callable[[float, float], float]\n",
        "ColorFunc3: TypeAlias = Callable[[int, float, float], float]\n",
        "ColorFunc4: TypeAlias = Callable[[int, int, float, float], float]\n",
        "ColorFunc: TypeAlias = ColorFunc2 | ColorFunc3 | ColorFunc4\n",
        "ColorSpec: TypeAlias = mplt.ColorType | ColorFunc\n",
        "\n",
        "\n",
        "def chaos_game(\n",
        "    *map: AffineMapLike | Sequence[AffineMapLike],\n",
        "    iterations: int = 50000,\n",
        "    warmups: int = 100,\n",
        "    x0: float = 1.0,\n",
        "    y0: float = 1.0,\n",
        "    s: float = 0.1,\n",
        "    c: ColorFunc | str = '#1f77b4',\n",
        "    cmap: matplotlib.colors.Colormap | str = 'viridis',\n",
        ") -> None:\n",
        "    \"\"\"カオスゲーム（chaos game）を実行する。\n",
        "\n",
        "    Args:\n",
        "        map (AffineMapLike or sequence of AffineMapLike): 変換行列と重み。\n",
        "        iterations (int): 反復回数。\n",
        "        warmups (int): 記録を始める前の慣らし反復の回数\n",
        "        x0 (float): 初期位置のx座標。\n",
        "        y0 (float): 初期位置のy座標。\n",
        "        s (float): マーカーポイントの大きさ。\n",
        "        c (ColorFunc or str): マーカーポイントの色。\n",
        "        cmap (Colormap or str): カラーマップ。\n",
        "    \"\"\"\n",
        "    weights = []\n",
        "    maps = []\n",
        "    for i, x in enumerate(map):\n",
        "        if not _parse_maps(x, weights, maps):\n",
        "            if _is_sequence_but_not_ndarray(x):\n",
        "                for j, y in enumerate(x):\n",
        "                    if not _parse_maps(y, weights, maps):\n",
        "                        msg = f'第{i + 1}引数のインデックス{j}の要素が不適: {y}'\n",
        "                        raise ValueError(msg)\n",
        "            else:\n",
        "                msg = f'第{i + 1}引数が不適: {x}'\n",
        "                raise ValueError(msg)\n",
        "    for w in weights:\n",
        "        if not np.isfinite(w):\n",
        "            msg = f'有限でない重みが与えられた: {w}'\n",
        "            raise ValueError(msg)\n",
        "        if w <= 0:\n",
        "            msg = f'0以下の重みが与えられた: {w}'\n",
        "            raise ValueError(msg)\n",
        "    if len(weights) == 0:\n",
        "        raise ValueError('1つ以上の写像が必要')\n",
        "    probs = weights / np.array(weights).sum()\n",
        "\n",
        "    xpoints = []\n",
        "    ypoints = []\n",
        "\n",
        "    r = Vec(x0, y0)\n",
        "\n",
        "    samples = np.random.choice(range(len(probs)), size=warmups, p=probs)\n",
        "    for i in samples:\n",
        "        r = maps[i] @ r\n",
        "\n",
        "    samples = np.random.choice(range(len(probs)), size=iterations, p=probs)\n",
        "    for i in samples:\n",
        "        r = maps[i] @ r\n",
        "        xpoints.append(r[0])\n",
        "        ypoints.append(r[1])\n",
        "\n",
        "    if callable(c):\n",
        "        sig = inspect.signature(c)\n",
        "        if len(sig.parameters) == 2:\n",
        "            cpoints = np.vectorize(c)(xpoints, ypoints)\n",
        "        elif len(sig.parameters) == 3:\n",
        "            cpoints = np.vectorize(c)(samples, xpoints, ypoints)\n",
        "        elif len(sig.parameters) == 4:\n",
        "            cpoints = np.vectorize(c)(\n",
        "                range(1, iterations + 1), samples, xpoints, ypoints\n",
        "            )\n",
        "        else:\n",
        "            cpoints=None\n",
        "            cmap = None\n",
        "    else:\n",
        "        cpoints = c\n",
        "        cmap = None\n",
        "\n",
        "    fig, ax = plt.subplots()\n",
        "    ax.scatter(xpoints, ypoints, s=s, c=cpoints, cmap=cmap, edgecolors='none')\n",
        "    ax.set_aspect('equal')\n",
        "    plt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "`R(theta)`: 回転行列\n",
        "\n",
        "`S(sx, sy)`: 拡大縮小行列\n",
        "\n",
        "`T(tx, ty)`: 平行移動行列\n",
        "\n",
        "`Sh(shx, shy)`: 剪断行列\n",
        "\n",
        "`A(a11, a12, a13, a21, a22, a23)`: 一般のアフィン行列\n",
        "\n",
        "`pi`: 円周率, `e`: 自然対数の底, `phi`: 黄金比\n",
        "\n",
        "`sqrt(x)`: `x`の平方根, `sin(x)`: `x`の正弦, `cos(x)`: `x`の余弦\n",
        "\n",
        "`@`: 行列積, `inv(m)`: 行列`m`の逆行列, `m.T`: 行列`m`の転置"
      ],
      "metadata": {
        "id": "paxmhA7qvyjI"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "m1 = R(pi / 6) @ S(1 / sqrt(3), -1 / sqrt(3))\n",
        "m2 = T(1 / 2, sqrt(3) / 6) @ R(-pi / 6) @ S(1 / sqrt(3), -1 / sqrt(3))\n",
        "display(m1, m2)\n",
        "chaos_game(m1, m2)"
      ],
      "metadata": {
        "id": "rBTml2QSGYUP"
      },
      "execution_count": null,
      "outputs": []
    }
  ]
}