シェルピンスキーのギャスケット

家のデスクトップで久しぶりにRStudioを立ち上げたら、
ワークスペースに面白いものが残っていたので改めて。
（RのコードはどこかのWebサイトを参考にしたもの）

gasket <- function(n=100000) {
 t <- proc.time()
 m <- matrix(c(0.5, 0, 0, 0.5) , 2)
 lst <- matrix(0, n, 2) # matrixを用意
 lst[1,] <- c(100, 0)   # 初期値

 prb <- sample(1:3, n,rep=TRUE)

 for(i in 2:n){
   lst[i,] <- switch(prb[i], # switchで振り分け
                     m %*% lst[i-1,],
                     m %*% lst[i-1,] + c(100,0),
                     m %*% lst[i-1,] + c(50, 30))
 }
 plot(lst, pch=".", axes=F, asp=sqrt(3)*1.5)
 timeval <- proc.time() - t
 print(timeval)
}

## 10年以上前のPC（多分Duronだった気がする） 
# > gasket()
#  user  system elapsed
#  0.94    1.07    2.29
## 今年買ったノートPC（Core i7 10th）
# > gasket()
#   user  system elapsed 
#   0.27    0.50    0.79  

# %matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
import time
import statistics


proc_time_l = []
for _ in range(100):
#856 ms ± 8.06 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
    s = time.perf_counter()
    n=100000
    m = np.array(
       [[0.5, 0],
       [0, 0.5]]
    )
    lst = np.zeros((n, 2))
    lst[0] = [100, 0]
    prb = np.random.randint(3, size=n)

    for i in range(1, n):
       if prb[i] == 0:
           lst[i] = np.dot(m, lst[i-1])
       elif prb[i] == 1:
           lst[i] = np.dot(m, lst[i-1]) + [100, 0]
       else:
           lst[i] = np.dot(m, lst[i-1]) + [50, 30]

    fig = plt.figure(figsize=(13, 10))
    plt.plot(lst[:, 0], lst[:, 1], 'k.')
    plt.axis('off')
    proc_time_l.append(time.perf_counter()-s)
print(statistics.mean(proc_time_l))
print(statistics.stdev(proc_time_l))

0.49452882845999396
0.014653359037954079

ちょっとだけ改良。

%%timeit
import numpy as np
import matplotlib.pyplot as plt


n=100000
m = np.array(
   [[0.5, 0],
   [0, 0.5]]
)
out = np.empty((n, 2))
a = np.array([[0, 0], [100, 0], [50, 30]])
r = np.random.randint(3, size=n-1)
out[0] = [100, 0]
out[1:] = a[r]

for i in range(1, n):
    out[i] = out[i] + m.dot(out[i-1])

fig = plt.figure(figsize=(13, 10))
plt.plot(out[:, 0], out[:, 1], 'k.')
plt.axis('off')

1 loop, best of 3: 240 ms per loop

行列積の（線形じゃ無い）部分が再帰構造的処理になっているので、ループは無くせない（と思う）。
単純なループのみが残っているので、NumbaやCythonが効果的だろう。

import numpy as np
import matplotlib.pyplot as plt
from numba import jit


@jit(nopython=True, parallel=True, fastmath=True)
def f(out, n):
    m = np.array([[0.5, 0], [0, 0.5]])
    for i in range(1, n):
        out[i] = out[i] + m.dot(out[i-1])

%%timeit
n=100000
out = np.empty((n, 2))
a = np.array([[0, 0], [100, 0], [50, 30]])
r = np.random.randint(3, size=n-1)
out[0] = [100, 0]
out[1:] = a[r]
f(out, n)
fig = plt.figure(figsize=(13, 10))
plt.plot(out[:, 0], out[:, 1], 'k.')
plt.axis('off')

1 loop, best of 3: 139 ms per loop

%%cython
import cython
import numpy as np
import matplotlib.pyplot as plt
cimport numpy as np


@cython.ccall
@cython.locals(out=np.ndarray, m=np.ndarray, i=cython.uint, n=cython.uint)
def f(out, n):
    m = np.array([[0.5, 0], [0, 0.5]])
    for i in range(1, n):
        out[i] = out[i] + m.dot(out[i-1])

@cython.ccall
@cython.locals(out=np.ndarray, a=np.ndarray, r=np.ndarray, n=cython.uint)
def main():
    n=100000
    out = np.empty((n, 2))
    a = np.array([[0, 0], [100, 0], [50, 30]])
    r = np.random.randint(3, size=n-1)
    out[0] = [100, 0]
    out[1:] = a[r]
    f(out, n)
    fig = plt.figure(figsize=(13, 10))
    plt.plot(out[:, 0], out[:, 1], 'k.')
    plt.axis('off')

%%timeit
main()

1 loop, best of 3: 203 ms per loop

ほぼCの様な書き方にしたらもっと速くなるが、Pythonのメリットを失うので、スケーラビリティを考えれば、ループ部分の高速化だけであればNumbaがリーズナブルチョイス。