3戦略の進化ゲームのプロットをさまざまな方法で行う。
3戦略の頻度をそれぞれ$x_1, x_2, x_3$とし、$x_1 + x_2+ x_3 = 1$である。

Contents

（前編）

matplotlibでの三角図の作成
進化ゲームによる3戦略のダイナミクスをベクトル場で表す
進化ゲームによる3戦略のダイナミクスをヌルクラインで表す
固定点をプロットする

（後編）

セパラトリクスのプロットする
シミュレーション結果をヒートマップで重ねる

matplotlibでの三角図の作成

まずは、必要なライブラリのimportを行う

from scipy.integrate import solve_ivp
import scipy.integrate as integrate
import numpy as np
import matplotlib.pyplot as plt
import math
import sympy as sym
import pandas as pd

基本となる三角図を作る。
通常のプロットではx-y軸を表す四角い枠ができる。

fig = plt.figure()
ax1 = fig.add_subplot(111)

この枠線の中に三角形を書くため、四角形の枠は以下のコードで消す。

plt.tick_params(labelbottom=False, labelleft=False,
                labelright=False, labeltop=False)
plt.tick_params(bottom=False, left=False, right=False, top=False)
plt.gca().spines['bottom'].set_visible(False)
plt.gca().spines['left'].set_visible(False)
plt.gca().spines['right'].set_visible(False)
plt.gca().spines['top'].set_visible(False)

三角形を書くコードは以下のとおり。

translation = 0.5
h = np.sqrt(3.0)*0.5
ax1.plot([1.0-translation, 0.5-translation],
         [0.0, h], 'black', zorder=10, lw=3)
ax1.plot([0.0-translation, 0.5-translation],
         [0.0, h], 'black', zorder=10, lw=3)
ax1.plot([0.0-translation, 1.0-translation],
         [0.0, 0.0], 'black', zorder=10, lw=3)

translationで三角形の位置を調整する(translation=0.5で、横軸の値が0のところが三角形の中心となる）。

四角形の枠は以下のコードで消す。

plt.tick_params(labelbottom=False, labelleft=False,
                labelright=False, labeltop=False)
plt.tick_params(bottom=False, left=False, right=False, top=False)
plt.gca().spines['bottom'].set_visible(False)
plt.gca().spines['left'].set_visible(False)
plt.gca().spines['right'].set_visible(False)
plt.gca().spines['top'].set_visible(False)

各頂点の座標と名前を表示する。

tex1= r'$(0, \frac{\sqrt{3}}{2})$'
tex2=r'$(-0.5,0)$'
tex3=r'$(0.5,0)$'
plt.text(0.15, h+0.05, tex1, fontsize=15)
plt.text(-0.75,-0.1, tex2, fontsize=15)
plt.text(0.6, -0.1, tex3, fontsize=15)


vertex1, vertex2, vertex3 = r'$x_1$',r'$x_2$', r'$x_3$'
plt.text(-0.07, h+0.05, vertex3, fontsize=15)
plt.text(-0.58, 0.01, vertex1, fontsize=15)
plt.text(0.5, 0.01, vertex3, fontsize=15)

これまでをまとめると、

fig = plt.figure()
ax1 = fig.add_subplot(111)
plt.tick_params(labelbottom=False, labelleft=False,
                labelright=False, labeltop=False)
plt.tick_params(bottom=False, left=False, right=False, top=False)
plt.gca().spines['bottom'].set_visible(False)
plt.gca().spines['left'].set_visible(False)
plt.gca().spines['right'].set_visible(False)
plt.gca().spines['top'].set_visible(False)

translation = 0.5
h = np.sqrt(3.0)*0.5
ax1.plot([1.0-translation, 0.5-translation],
         [0.0, h], 'black', zorder=10, lw=3)
ax1.plot([0.0-translation, 0.5-translation],
         [0.0, h], 'black', zorder=10, lw=3)
ax1.plot([0.0-translation, 1.0-translation],
         [0.0, 0.0], 'black', zorder=10, lw=3)

tex1= r'$(0, \frac{\sqrt{3}}{2})$'
tex2=r'$(-0.5,0)$'
tex3=r'$(0.5,0)$'
plt.text(0.15, h+0.05, tex1, fontsize=15)
plt.text(-0.75,-0.1, tex2, fontsize=15)
plt.text(0.6, -0.1, tex3, fontsize=15)

vertex1, vertex2, vertex3 = r'$x_1$',r'$x_2$', r'$x_3$'
plt.text(-0.07, h+0.05, vertex3, fontsize=15)
plt.text(-0.58, 0.01, vertex1, fontsize=15)
plt.text(0.5, 0.01, vertex3, fontsize=15)
plt.show()

三角形の各頂点はある一つの戦略が集団を占めている状況を意味する。

進化ゲームによる3戦略のダイナミクスをベクトル場で表す

いよいよベクトル場を書くが、その前に、どのような変換をすれば3次元の値を２次元に落とし込めるか考えてみたい。
３次元の値を2次元にプロットするため、正三角形と次のようなベクトルを考えてみよう。

座標$(0,0)$から各頂点へ向かう3つのベクトル$\vec{e_1},\vec{e_2}, \vec{e_3}$は、二次元の座標に変換した時のそれぞれの戦略の頻度を表している。各ベクトルは、

\begin{align}
 \vec{e_1} &= (-0.5, 0)\\
 \vec{e_2} &= (0.5,0)\\
\vec{e_3}&= (0,\frac{\sqrt{3}}{2})
\end{align}

であり、三角図内の任意の位置ベクトル$\vec{a}$$(u, v)$はこれらのベクトルを使って以下のように表せる。

 \vec{a} = x_1 \vec{e_1} + x_2 \vec{e_2}+ x_3 \vec{e_3}

ここにベクトルの値を入れて計算すると、

 （u, v) = (-0.5 x_1+0.5x_2, \frac{\sqrt{3}}{2}x_3) \tag{1}

となり、二次元で3つの戦略の頻度を表すことができる。

上記の議論を参考に、ベクトル場をプロットする。３次元でのダイナミクスを計算し、その後、２次元に変換したものを三角図にプロットする方針で進める。
最初に、レプリケーターダイナミクスを定義する。

#利得行列の要素
a11,a12,a13= 1,0,0
a21,a22,a23=0,1,0
a31,a32,a33=0,0,1

def dx(x, y):
    z = 1-x-y
    f_A = a11*x+ a12*y+a13*z
    f_B = a21*x+ a22*y+a23*z
    f_C = a31*x+ a32*y+a33*z
    phi = x*f_A + y*f_B + z*f_C
    return x*(f_A - phi)

def dy(x, y):
    z = 1-x-y
    f_A = a11*x+ a12*y+a13*z
    f_B = a21*x+ a22*y+a23*z
    f_C = a31*x+ a32*y+a33*z
    phi = x*f_A + y*f_B + z*f_C
    return y*(f_B - phi)

ここからのコードでは、$x_1, x_2,x_3$の表記を変えて、$x,y,z$でそれぞれの頻度を表している。
zの頻度は1-x-yで表され、その時間微分は$-\dot{x}-\dot{y}$で計算できる。
ベクトル場は、quiver関数を使ってプロットする。
・quiver($ X, Y,\Delta X,\Delta V$)
$X, Y$がベクトルの開始点、$\Delta X,\Delta V$は変化量である。行列を使って指定した座標上の全体から複数の矢印をプロットできる。
ベクトルを計算するための３次元のグリッドを作成する。この際、仮定していない領域（$z=1-x-y<0$)はnoneを代入し除外する。

# 作図範囲を設定するパラメータ
boundary = 0.005

# 三次元座標のXの計算する範囲は[0, 1]であるが、二次元座標のUの場合は[-0.5, 0.5]であるため、プロットする範囲はUの範囲に設定する
xmax, xmin = 1, 0.01
ymax, ymin = 1, 0.01

v_X, v_Y = np.meshgrid(np.arange(xmin + 9*boundary, ymax-9*boundary, 0.05),
                       np.arange(xmin + 9*boundary, ymax - 9*boundary, 0.05))
for i in range(len(v_X)):
    for j in range(len(v_X)):
        if v_X[i][j] + v_Y[i][j] > 1.006:
            v_X[i][j] = None
            v_Y[i][j] = None

これでベクトル場用の三次元のメッシュ$v_X, v_Y, v_Z(=1-v_X - v_Y)$ができた。これを上記の変換式$(1)$を使って二次元のメッシュ$(v_U, v_V)$に変換する。

v_U = -0.5*v_X + 0.5*v_Y
v_V = (1 - v_X - v_Y)*h

微分方程式も同様に計算し、最後にquiverでプロットする。quiverの３,4番目の引数はベクトルの大きさで割っている。こうすることでベクトルの大きさを揃えている(R)。
図の透明度はalphaで指定する。

v_dX = dx(v_X, v_Y)
v_dY = dy(v_X, v_Y)
v_dU = -0.5*v_dX + 0.5*v_dY
v_dV = -v_dX - v_dY
plt.quiver(v_U, v_V, v_dU/np.sqrt(pow(v_dU, 2)+pow(v_dV, 2)),v_dV/np.sqrt(pow(v_dU, 2)+pow(v_dV, 2)), alpha=0.5)

これまでをまとめると

fig = plt.figure()
ax1 = fig.add_subplot(111)
plt.tick_params(labelbottom=False, labelleft=False,
                labelright=False, labeltop=False)
plt.tick_params(bottom=False, left=False, right=False, top=False)
plt.gca().spines['bottom'].set_visible(False)
plt.gca().spines['left'].set_visible(False)
plt.gca().spines['right'].set_visible(False)
plt.gca().spines['top'].set_visible(False)

translation = 0.5
h = np.sqrt(3.0)*0.5
ax1.plot([1.0-translation, 0.5-translation],
         [0.0, h], 'black', zorder=10, lw=3)
ax1.plot([0.0-translation, 0.5-translation],
         [0.0, h], 'black', zorder=10, lw=3)
ax1.plot([0.0-translation, 1.0-translation],
         [0.0, 0.0], 'black', zorder=10, lw=3)

tex1= r'$(0, \frac{\sqrt{3}}{2})$'
tex2=r'$(-0.5,0)$'
tex3=r'$(0.5,0)$'
plt.text(0.15, h+0.05, tex1, fontsize=15)
plt.text(-0.75,-0.1, tex2, fontsize=15)
plt.text(0.6, -0.1, tex3, fontsize=15)

vertex1, vertex2, vertex3 = r'$x_1$',r'$x_2$', r'$x_3$'
plt.text(-0.07, h+0.05, vertex3, fontsize=15)
plt.text(-0.58, 0.01, vertex1, fontsize=15)
plt.text(0.5, 0.01, vertex3, fontsize=15)

# 利得行列の成分
a11,a12,a13= 1,0,0
a21,a22,a23=0,1,0
a31,a32,a33=0,0,1

def dx(x, y):
    z = 1-x-y
    f_A = a11*x+ a12*y+a13*z
    f_B = a21*x+ a22*y+a23*z
    f_C = a31*x+ a32*y+a33*z
    phi = x*f_A + y*f_B + z*f_C
    return x*(f_A - phi)


def dy(x, y):
    z = 1-x-y
    f_A = a11*x+ a12*y+a13*z
    f_B = a21*x+ a22*y+a23*z
    f_C = a31*x+ a32*y+a33*z
    phi = x*f_A + y*f_B + z*f_C
    return y*(f_B - phi)

def dz(x, y):
    z = 1-x-y
    f_A = a11*x+ a12*y+a13*z
    f_B = a21*x+ a22*y+a23*z
    f_C = a31*x+ a32*y+a33*z
    phi = x*f_A + y*f_B + z*f_C
    return z*(f_C - phi)

# 作図範囲を設定するパラメータ
p = 2.0
boundary = 0.005

# 三次元座標のXの計算する範囲は[0, 1]であるが、二次元座標のUの場合は[-0.5, 0.5]であるため、プロットする範囲はUの範囲に設定する
xmax, xmin = 1, 0.01
ymax, ymin = 1, 0.01
umax, umin = 0.5, -0.5

v_X, v_Y = np.meshgrid(np.arange(xmin + 9*boundary, ymax-9*boundary, 0.05),
                       np.arange(xmin + 9*boundary, ymax - 9*boundary, 0.05))
for i in range(len(v_X)):
    for j in range(len(v_X)):
        if v_X[i][j] + v_Y[i][j] > 1.006:
            v_X[i][j] = None
            v_Y[i][j] = None

v_U = -0.5*v_X + 0.5*v_Y
v_V = (1 - v_X - v_Y)*h

v_dX = dx(v_X, v_Y)
v_dY = dy(v_X, v_Y)
v_dU = -0.5*v_dX + 0.5*v_dY
v_dV = (-v_dX - v_dY)
a = plt.quiver(v_U, v_V, v_dU/np.sqrt(pow(v_dU, 2)+pow(v_dV, 2)),v_dV/np.sqrt(pow(v_dU, 2)+pow(v_dV, 2)), alpha=0.5)

plt.show()

進化ゲームによる3戦略のダイナミクスをヌルクラインで表す

基本的なやり方はベクトルと同じで各座標を表すmesh gridを用意する。これをquiverのかわりにcontour関数を使うことで、それぞれの変数の時間微分が0となっている線を描く。
・contour([X, Y,] Z, [levels], **kwargs)
X,Yは値Zの座標(x,y)と対応している。levels=0とすると、1つの曲線が引かれ、これが時間微分=0と対応するようだ（R)。

U = -0.5 * X + 0.5*Y
V = (1-X-Y)*h
dX = dx(X, Y)
dY = dy(X, Y)
dZ = -dx(X, Y)-dy(X, Y)
dU = -0.5*dX + 0.5*dY
dV = dZ*h

plt.legend(loc='upper left', bbox_to_anchor=(1, 1))

plt.contour(U, V,  dY, levels=[0], colors="red",
            linewidths=2, )
plt.contour(U, V, dX, levels=[0], colors="Blue",
            linewidths=2, )
plt.contour(U, V,  dZ, levels=[0], colors="green", linewidths=2)
# 凡例用(quiverではおそらく凡例が表示されないため、かわりにplotで作成する。plotは表示範囲外にするために、xlim,ylimで表示範囲を制限する。)
ax1.plot([10, 10], 'blue', lw=1, label="x ")
ax1.plot([11, 11], 'red', lw=1, label="y ")
ax1.plot([12, 12], 'green', lw=1, label="z ")
ax1.legend(loc='upper left', bbox_to_anchor=(1, 1))
plt.xlim(-0.7, 0.7)
plt.ylim(-0.1, 1)

固定点をプロットする

次のような利得行列を考える。

a11,a12,a13= 1, 0, 0.6
a21,a22,a23= 1, 1, 0.6
a31,a32,a33=0.3, 0.3, 2

このとき、図はこのようになる。

このときの固定点を安定であれば黒色に、不安定であれば白色にプロットしたい。
方法としては、sympyを使って各固定点の安定性を計算したあと、色を塗り分けるため安定性により条件分岐し、scatterでplotする。

 変数、定数を宣言
x = sym.Symbol('x')
y = sym.Symbol('y')

# 微分方程式の定義
z = 1-x-y
f_A = a11*x+ a12*y+a13*z
f_B = a21*x+ a22*y+a23*z
f_C = a31*x+ a32*y+a33*z
phi = x*f_A + y*f_B + (1-x-y)*f_C

dx = x*(f_A - phi)
dy = y*(f_B - phi)

# 固定点の座標を求める
eq1 = sym.Eq(dx, 0)
eq2 = sym.Eq(dy, 0)
equivs = sym.solve([eq1, eq2], [x, y])

# ヤコビ行列の各要素を求める
j11 = sym.diff(dx, x)
j12 = sym.diff(dx, y)
j21 = sym.diff(dy, x)
j22 = sym.diff(dy, y)

# 要素を並べて行列にする
J = sym.Matrix(([j11, j12], [j21, j22]))
u = sym.Matrix((x, y))


# それぞれの平衡点についての固有値を求める
for i in range(len(equivs)):
    xy_equiv1 = list(zip(u, equivs[i]))
    ev1 = J.subs(xy_equiv1).eigenvals()

    for lam, dup in ev1.items():
        # ここで、固有値が重複している場合にも対応できるようにする。
        tr = tr + lam*dup
        det = det*lam**dup

    cordinate_u = -0.5*equivs[i][0]+0.5*equivs[i][1]
    cordinate_v = (1 - equivs[i][0] - equivs[i][1])*h

# tr, detが複素数であった場合にも対応
    if cordinate_v >= -0.01:
        if sym.re(det) > 0 and sym.re(tr) < 0 and sym.im(cordinate_u) == 0 and sym.im(cordinate_v) == 0:
            ax1.scatter(
                cordinate_u, cordinate_v,
                zorder=10,
                s=100,
                c='black',
                marker='o',
                linewidths=1,
                edgecolors='black')
        elif sym.im(cordinate_u) == 0 and sym.im(cordinate_v) == 0 and cordinate_v < h+0.1 and cordinate_v >= 0.0:
            ax1.scatter(
                cordinate_u, cordinate_v,
                zorder=10,
                s=100,
                c='white',
                marker='o',
                linewidths=1,
                edgecolors='black')

これまでのコードを以下にまとめた。

fig = plt.figure()
ax1 = fig.add_subplot(111)
plt.tick_params(labelbottom=False, labelleft=False,
                labelright=False, labeltop=False)
plt.tick_params(bottom=False, left=False, right=False, top=False)
plt.gca().spines['bottom'].set_visible(False)
plt.gca().spines['left'].set_visible(False)
plt.gca().spines['right'].set_visible(False)
plt.gca().spines['top'].set_visible(False)

translation = 0.5
h = np.sqrt(3.0)*0.5
ax1.plot([1.0-translation, 0.5-translation],
         [0.0, h], 'black', zorder=10, lw=3)
ax1.plot([0.0-translation, 0.5-translation],
         [0.0, h], 'black', zorder=10, lw=3)
ax1.plot([0.0-translation, 1.0-translation],
         [0.0, 0.0], 'black', zorder=10, lw=3)

tex1= r'$(0, \frac{\sqrt{3}}{2})$'
tex2=r'$(-0.5,0)$'
tex3=r'$(0.5,0)$'
plt.text(0.15, h+0.05, tex1, fontsize=15)
plt.text(-0.75,-0.1, tex2, fontsize=15)
plt.text(0.6, -0.1, tex3, fontsize=15)

vertex1, vertex2, vertex3 = r'$x_1$',r'$x_2$', r'$x_3$'
plt.text(-0.07, h+0.05, vertex3, fontsize=15)
plt.text(-0.58, 0.01, vertex1, fontsize=15)
plt.text(0.5, 0.01, vertex2, fontsize=15)

a11,a12,a13= 1, 0, 0.6
a21,a22,a23= 1, 1.1, 0.6
a31,a32,a33=0.3, 0.3, 2

def dx(x, y):
    z = 1-x-y
    f_A = a11*x+ a12*y+a13*z
    f_B = a21*x+ a22*y+a23*z
    f_C = a31*x+ a32*y+a33*z
    phi = x*f_A + y*f_B + z*f_C
    return x*(f_A - phi)


def dy(x, y):
    z = 1-x-y
    f_A = a11*x+ a12*y+a13*z
    f_B = a21*x+ a22*y+a23*z
    f_C = a31*x+ a32*y+a33*z
    phi = x*f_A + y*f_B + z*f_C
    return y*(f_B - phi)

def dz(x, y):
    z = 1-x-y
    f_A = a11*x+ a12*y+a13*z
    f_B = a21*x+ a22*y+a23*z
    f_C = a31*x+ a32*y+a33*z
    phi = x*f_A + y*f_B + z*f_C
    return z*(f_C - phi)

# 作図範囲を設定するパラメータ
p = 2.0
boundary = 0.005

# 三次元座標のXの計算する範囲は[0, 1]であるが、二次元座標のUの場合は[-0.5, 0.5]であるため、プロットする範囲はUの範囲に設定する
xmax, xmin = 1, 0.01
ymax, ymin = 1, 0.01
umax, umin = 0.5, -0.5

v_X, v_Y = np.meshgrid(np.arange(xmin + 9*boundary, ymax-9*boundary, 0.05),
                       np.arange(xmin + 9*boundary, ymax - 9*boundary, 0.05))
for i in range(len(v_X)):
    for j in range(len(v_X)):
        if v_X[i][j] + v_Y[i][j] > 1.006:
            v_X[i][j] = None
            v_Y[i][j] = None

v_U = -0.5*v_X + 0.5*v_Y
v_V = (1 - v_X - v_Y)*h

v_dX = dx(v_X, v_Y)
v_dY = dy(v_X, v_Y)
v_dU = -0.5*v_dX + 0.5*v_dY
v_dV = (-v_dX - v_dY)
a = plt.quiver(v_U, v_V, v_dU/np.sqrt(pow(v_dU, 2)+pow(v_dV, 2)),v_dV/np.sqrt(pow(v_dU, 2)+pow(v_dV, 2)), alpha=0.5)


X, Y = np.meshgrid(np.arange(xmin - boundary, xmax + boundary, 0.005),
                   np.arange(ymin - boundary, ymax + boundary, 0.005))
for i in range(len(X)):
    for j in range(len(X)):
        if X[i][j] + Y[i][j] > 1.006:
            X[i][j] = None
            Y[i][j] = None


U = -0.5 * X + 0.5*Y
V = (1-X-Y)*h
dX = dx(X, Y)
dY = dy(X, Y)
# dZ = dz(X, Y)
dZ = -dx(X, Y)-dy(X, Y)
dU = -0.5*dX + 0.5*dY
dV = dZ*h

plt.legend(loc='upper left', bbox_to_anchor=(1, 1))

plt.contour(U, V,  dY, levels=[0], colors="red",
            linewidths=2, )
plt.contour(U, V, dX, levels=[0], colors="Blue",
            linewidths=2, )
plt.contour(U, V,  dZ, levels=[0], colors="green", linewidths=2)
# 凡例用(quiverではおそらく判例が表示されないため、かわりにplotで作成する。plotは表示範囲外にするために、xlim,ylimで表示範囲を制限する)
ax1.plot([10, 10], 'blue', lw=1, label="x ")
ax1.plot([11, 11], 'red', lw=1, label="y ")
ax1.plot([12, 12], 'green', lw=1, label="z ")
ax1.legend(loc='upper left', bbox_to_anchor=(1, 1))
plt.xlim(-0.7, 0.7)
plt.ylim(-0.1, 1)


# 固定点のプロット
# 固定点の座標を求める、固定点の安定性を調べる、安定性によって点の色を変えるという流れで進む。

# 変数、定数を宣言
x = sym.Symbol('x')
y = sym.Symbol('y')

# 微分方程式の定義
z = 1-x-y
f_A = a11*x+ a12*y+a13*z
f_B = a21*x+ a22*y+a23*z
f_C = a31*x+ a32*y+a33*z
phi = x*f_A + y*f_B + (1-x-y)*f_C

dx = x*(f_A - phi)
dy = y*(f_B - phi)

# 固定点の座標を求める
eq1 = sym.Eq(dx, 0)
eq2 = sym.Eq(dy, 0)
equivs = sym.solve([eq1, eq2], [x, y])

# ヤコビ行列の各要素を求める
j11 = sym.diff(dx, x)
j12 = sym.diff(dx, y)
j21 = sym.diff(dy, x)
j22 = sym.diff(dy, y)

# 要素を並べて行列にする
J = sym.Matrix(([j11, j12], [j21, j22]))
u = sym.Matrix((x, y))


# それぞれの平衡点についての固有値を求める
for i in range(len(equivs)):
    xy_equiv1 = list(zip(u, equivs[i]))
    ev1 = J.subs(xy_equiv1).eigenvals()

    for lam, dup in ev1.items():
        # ここで、固有値が重複している場合にも対応できるようにする。
        tr = tr + lam*dup
        det = det*lam**dup

    cordinate_u = -0.5*equivs[i][0]+0.5*equivs[i][1]
    cordinate_v = (1 - equivs[i][0] - equivs[i][1])*h

# tr, detが複素数であった場合にも対応
    if cordinate_v >= -0.01:
        if sym.re(det) > 0 and sym.re(tr) < 0 and sym.im(cordinate_u) == 0 and sym.im(cordinate_v) == 0:
            ax1.scatter(
                cordinate_u, cordinate_v,
                zorder=10,
                s=100,
                c='black',
                marker='o',
                linewidths=1,
                edgecolors='black')
        elif sym.im(cordinate_u) == 0 and sym.im(cordinate_v) == 0 and cordinate_v < h+0.1 and cordinate_v >= 0.0:
            ax1.scatter(
                cordinate_u, cordinate_v,
                zorder=10,
                s=100,
                c='white',
                marker='o',
                linewidths=1,
                edgecolors='black')
plt.show()

引用、謝辞、その他備考

上記のコードは、昔にさまざまなサイトを参考にして書きました。ここでそれらの記事の投稿者様に感謝を申し上げます。覚えている範囲で参考にしたのはからっぽの書庫とQiitaの記事、そうだ！研究しようです(敬称略）。
他のサイトの投稿者様については失念してしまいましたが、もし似たような書き方をしているサイトがございましたら、コメントをいただけると幸いです。私が参考にしたと思われる場合、引用、謝辞に追加したいと思います。

後編は、気が乗ったら書きます。