昨天我看到了關于牛頓分形的新 3Blue1Brown 視頻，我真的被他對分形的現場表現所迷住了。（這是任何感興趣的人的視頻鏈接，它在 13:40：https : //www.youtube.com/watch?v= -RdOwhmqP5s ）
我想自己嘗試一下并嘗試用 python 撰寫它（我認為他也使用 python）。
我花了幾個小時試圖改進我的幼稚實作，但到了我不知道如何讓它更快的地步。
代碼如下所示：
import os
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.gridspec import GridSpec
from time import time
def print_fractal(state):
fig = plt.figure(figsize=(8, 8))
gs = GridSpec(1, 1)
axs = [fig.add_subplot(gs[0, 0])]
fig.tight_layout(pad=5)
axs[0].matshow(state)
axs[0].set_xticks([])
axs[0].set_yticks([])
plt.show()
plt.close()
def get_function_value(z):
return z**5 z**2 - z 1
def get_function_derivative_value(z):
return 5*z**4 2*z - 1
def check_distance(state, roots):
roots2 = np.zeros((roots.shape[0], state.shape[0], state.shape[1]), dtype=complex)
for r in range(roots.shape[0]):
roots2[r] = np.full((state.shape[0], state.shape[1]), roots[r])
dist_2 = np.abs((roots2 - state))
original_state = np.argmin(dist_2, axis=0) 1
return original_state
def static():
time_start = time()
s = 4
c = [0, 0]
n = 800
polynomial = [1, 0, 0, 1, -1, 1]
roots = np.roots(polynomial)
state = np.transpose((np.linspace(c[0] - s/2, c[0] s/2, n)[:, None] 1j*np.linspace(c[1] - s/2, c[1] s/2, n)))
n_steps = 15
time_setup = time()
for _ in range(n_steps):
state -= (get_function_value(state) / get_function_derivative_value(state))
time_evolution = time()
original_state = check_distance(state, roots)
time_check = time()
print_fractal(original_state)
print("{0:<40}".format("Time to setup the initial configuration:"), "{:20.3f}".format(time_setup - time_start))
print("{0:<40}".format("Time to evolve the state:"), "{:20.3f}".format(time_evolution - time_setup))
print("{0:<40}".format("Time to check the closest roots:"), "{:20.3f}".format(time_check - time_evolution))
平均輸出如下所示：
設定初始配置的時間：0.004
狀態進化時間：0.796
檢查最近根的時間：0.094
It's clear that it's the evolution part that bottlenecks the process. It's not "slow", but I think it's not enough to render something live like in the video. I already did what I could by using numpy vectors and avoiding loops but I guess it's not enough. What other tricks could be applied here?
Note: I tried using numpy.polynomials.Polynomial class to evaluate the function, but it was slower than this version.
uj5u.com熱心網友回復：
for _ in range(n_steps):
state -= (get_function_value(state) / get_function_derivative_value(state))
如果你有足夠的記憶體，你可以嘗試向量化這個回圈并使用矩陣計算存盤每個迭代步驟。
uj5u.com熱心網友回復：
- 通過使用單個復數 (
np.complex64) 精度，我得到了改進（快了約 40% ）。
(...)
state = np.transpose((np.linspace(c[0] - s/2, c[0] s/2, n)[:, None]
1j*np.linspace(c[1] - s/2, c[1] s/2, n)))
state = state.astype(np.complex64)
(...)
- 3Blue1Brown 在描述中添加了這個鏈接：https ://codepen.io/mherreshoff/pen/RwZPazd 你可以看看它是如何完成的（旁注：這支筆的作者也使用了單精度）
轉載請註明出處，本文鏈接：https://www.uj5u.com/houduan/325915.html
標籤：python numpy performance optimization fractals
上一篇：概率變異資料幀