家のデスクトップで久しぶりに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がリーズナブルチョイス。