np.where

np.where은 특정 조건에 따라 ndarray에 연산을 행하여 반환합니다.

numpy.where
Return elements chosen from x or y depending on condition.

이는 두 가지의 기능으로 구분할 수 있습니다.

첫 번째로, 어떤 ndarray에 조건 연산을 취한 boolean 배열을 인자로 넣음으로써 이를 만족하는 배열의 원소들의 index를 반환합니다. 가령,

np.where(np.arange(100) % 2 == 0)

일 경우, np.arange(100)의 원소 중 짝수인 것들의 index를 반환하므로 반환값은 ndarray([0, 2, ..., 96, 98]) 이 됩니다.

두 번째로, 이렇게 조건에 따라 참인 것과 거짓인 원소에 대해 각기 다른 연산을 시행하여 그 결과값을 반환합니다. 가령,

arr = np.arange(100)
np.where(arr % 2 == 0, arr, -arr)

일 경우, np.arange(100)의 원소 중 짝수인 것들은 그대로, 홀수인 것들은 그것에 음수 기호를 붙인 수를 반환하므로 반환값은 ndarray([0, -1, 2, -3, ..., 98, -99])가 됩니다.

NumPy를 이용하여 Mandelbrot 그림 그리기

Mandelbrot 집합은 다음의 점화식으로 정의된 수열이 무한대로 발산하지 않게 하는 (진동하는 것은 가능합니다) 복소수 $c$의 집합으로 정의됩니다.

이에 대해서는 자기 유사성(프랙탈) 등 특이적인 성질을 가지고 있으나, 정작 복소수 $c$의 범위는 알 수 없으며 계산을 이용해 구하여야 합니다. 수열이 무한대로 발산하지 않는다라는 조건 때문에 정확히 집합을 그리기 위해서는 무한 번의 연산을 수행해야 할 것이나, 실제로는 그렇게 하는 것이 불가능하기 때문에 적당한 반복 회수에 걸쳐 발산 여부를 검증합니다. 그 방법 중 하나는 $z_n$의 modulus가 2보다 크다면 $(z_n)$은 발산한다는 성질을 이용해 반복 이후 $\vert z_n\vert>2$라면 발산함으로 가정합니다.

따라서 그림은 발산 여부 만으로 이루어진 점으로 구성되므로 흑백이 이상적이나, 보통 인터넷에서 볼 수 있는 그림은 발산 속도를 기준으로 색을 나누어 채색한 것입니다. 여기에서는 발산 여부만을 가지고 그림을 그려보도록 하겠습니다.

Mandelbrot 그림의 알려진 상하좌우의 경계는 [-2, 1]×[-1.5, 1.5]입니다. 먼저 연산과 시각화에 필요한 모듈부터 로드하도록 하겠습니다.

import numpy as np
import matplotlib.pyplot as plt

영역 [-2, 1]×[-1.5, 1.5]을 각각 가로/세로로 1000분할하여 미소 영역(점)마다 발산 여부를 계산하도록 하겠습니다. np.meshgrid를 이용하여 주어진 x좌표, y좌표의 Cartesian 곱을 구하겠습니다.

size = 1000

x = np.linspace(-2, 1, size+1)
y = np.linspace(-1.5, 1.5, size+1)
xx, yy = np.meshgrid(x, y)

모든 각 미소 영역(점)을 $c$로 하여 발산 여부를 결정해야 하므로, 구한 곱집합을 복소평면의 원소로 변환한 후 2차원 행렬화하겠습니다.

c = xx + yy * 1j # c는 1000x1000 행렬

또 각 $c$에 대해 $z$의 값도 계산해야 하므로 초기 z를 c와 같은 크기인 영행렬 (여기서 행렬의 원소는 복소수입니다)로 만들겠습니다.

z = np.zeros_like(c, dtype=complex) 

이제 각 $c$에 대한 발산 여부를 검증해야 합니다. 여기에서는 주어진 반복 횟수만큼 $z^2+c\mapsto z$를 반복하고, 최종적으로 $z$의 modulus가 2 이하인지 검증하도록 하겠습니다. 이는 최적화된 알고리즘은 아닙니다. 이보다 좋은 방법은 각각의 $c$에 대해 $z^2+c\mapsto z$의 modulus가 2보다 커지는 순간 중지하는 것이 계산량을 줄이는 방법입니다. 그러나 Python은 elementwise한 for 연산에 취약하고, 행렬 연산을 이용하는 것이 NumPy의 특장점인 만큼 조금 성능은 떨어지더라도 가독성이 좋은 방법을 사용하겠습니다.

상기 np.where 방법을 이용하면 $z$에 $z^2+c$를 대입하는 방법에서 약간의 overhead를 줄일 수 있습니다. $z$가 발산한다면 그 수의 절대값이 매우 커질 것이므로 발산하기 시작하는 시점(modulus가 2 이상인 시점)에서부터 계산을 하지 않고 그대로 $z$를 대입하는 방법을 고안할 수 있습니다. 따라서 다음과 같이 반복문을 구상할 수 있습니다.

for _ in range(100): # 100회 반복
    z = np.where(np.abs(z) < 2, z**2+c, z)

코드를 해석하면 다음과 같습니다.

마지막으로 반복을 마친 후 z에 대해 발산 여부(np.abs(z) < 2)를 가지고 boolean 행렬로 만든 후, plt.imshow를 이용해 이를 흑백 이미지로 출력하면 됩니다.

plt.imshow(np.abs(z) < 2, extent=[-2, 1, -1.5, 1.5])
plt.gray()
plt.show()