跳至內容

File:Exponential Collatz Fractal.jpg

頁面內容不支援其他語言。
這個檔案來自維基共享資源
維基百科,自由的百科全書

原始檔案 (30,720 × 17,280 像素,檔案大小:41.13 MB,MIME 類型:image/jpeg


摘要

警告 部分瀏覽器在瀏覽此圖片的完整大小時可能會遇到困難:該圖片中有數量巨大的像素點,可能無法完全載入或者導致您的瀏覽器停止回應。 交互式大图查看器
描述
English: A Collatz fractal for the interpolating function . The center of the image is and the real part goes from to .
日期
來源 自己的作品
作者 Hugo Spinelli
其他版本
new file This image is a JPEG version of the original PNG image at File: Exponential Collatz Fractal.png.

Generally, this JPEG version should be used when displaying the file from Commons, in order to reduce the file size of thumbnail images. However, any edits to the image should be based on the original PNG version in order to prevent generation loss, and both versions should be updated. Do not make edits based on this version.  See here for more information.

Source code
InfoField
Python code
import enum
import time

import numba as nb
import numpy as np
import matplotlib
from PIL import Image, PyAccess


# Amount of times to print the total progress
PROGRESS_STEPS: int = 20

# Number of pixels (width*height) and aspect ratio (width/height)
RESOLUTION: int = 1920*1080
ASPECT_RATIO: float = 1920/1080

# Value of the center pixel
CENTER: complex = 0 + 0j  # For testing: -4.6875 + 2.63671875j

# Value range of the real part (width of the horizontal axis)
RE_RANGE: float = 10  # For testing: 10/16

# Show grid lines for integer real and imaginary parts
SHOW_GRID: bool = False
GRID_COLOR: tuple[int, int, int] = (125, 125, 125)

# Matplotlib named colormap
COLORMAP_NAME: str = 'inferno'

# Color of the interior of the fractal (convergent points)
INTERIOR_COLOR: tuple[int, int, int] = (0, 0, 60)
# Color for large divergence counts
ERROR_COLOR: tuple[int, int, int] = (0, 0, 60)


# Plot range of the axes
X_MIN = CENTER.real - RE_RANGE/2  # min Re(z)
X_MAX = CENTER.real + RE_RANGE/2  # max Re(z)
Y_MIN = CENTER.imag - RE_RANGE/(2*ASPECT_RATIO)  # min Im(z)
Y_MAX = CENTER.imag + RE_RANGE/(2*ASPECT_RATIO)  # max Im(z)

x_range = X_MAX - X_MIN
y_range = Y_MAX - Y_MIN
pixels_per_unit = np.sqrt(RESOLUTION/(x_range*y_range))

# Width and height of the image in pixels
WIDTH = round(pixels_per_unit*x_range)
HEIGHT = round(pixels_per_unit*y_range)


# Maximum iterations for the divergence test
MAX_ITER: int = 10**4  # recommended >= 10**3
# Minimum consecutive abs(r) decreases to declare linear divergence
MIN_R_DROPS: int = 4  # recommended >= 2
# Minimum iterations to start checking for slow drift (unknown divergence)
MIN_ITER_SLOW: int = 200  # recommended >= 100


# Max value of Re(z) and Im(z) such that the recursion doesn't overflow
CUTOFF_RE = 7.564545572282618e+153
CUTOFF_IM = 112.10398935569289


# Precompute the colormap
CMAP_LEN: int = 2000
cmap_mpl = matplotlib.colormaps[COLORMAP_NAME]
# Start away from 0 (discard black values for the 'inferno' colormap)
# Matplotlib's colormaps have 256 discrete color points
n_cmap = round(256*0.85)
CMAP = [cmap_mpl(k/256) for k in range(256 - n_cmap, 256)]
# Interpolate
x = np.linspace(0, 1, num=CMAP_LEN)
xp = np.linspace(0, 1, num=n_cmap)
c0, c1, c2 = tuple(np.interp(x, xp, [c[k] for c in CMAP]) for k in range(3))
CMAP = []
for x0, x1, x2 in zip(c0, c1, c2):
    CMAP.append(tuple(round(255*x) for x in (x0, x1, x2)))


class DivType(enum.Enum):
    """Divergence type."""

    MAX_ITER = 0  # Maximum iterations reached
    SLOW = 1  # Detected slow growth (maximum iterations will be reached)
    CYCLE = 2  # Cycled back to the same value after 8 iterations
    LINEAR = 3  # Detected linear divergence
    CUTOFF_RE = 4  # Diverged by exceeding the real part cutoff
    CUTOFF_IM = 5  # Diverged by exceeding the imaginary part cutoff


@nb.jit(nb.float64(nb.float64, nb.int64), nopython=True)
def smooth(x, k=1):
    """Recursive exponential smoothing function."""

    y = np.expm1(np.pi*x)/np.expm1(np.pi)
    if k <= 1:
        return y
    return smooth(y, np.fmin(6, k - 1))


@nb.jit(nb.float64(nb.float64), nopython=True)
def get_delta_im(x):
    """Get the fractional part of the smoothed divergence count for
    imaginary part blow-up."""

    nu = np.log(np.abs(x)/CUTOFF_IM)/(np.pi*CUTOFF_IM - np.log(CUTOFF_IM))
    nu = np.fmax(0, np.fmin(nu, 1))
    return smooth(1 - nu, 2)


@nb.jit(nb.float64(nb.float64, nb.float64), nopython=True)
def get_delta_re(x, e):
    """Get the fractional part of the smoothed divergence count for
    real part blow-up."""

    nu = np.log(np.abs(x)/CUTOFF_RE)/np.log1p(e)
    nu = np.fmax(0, np.fmin(nu, 1))
    return 1 - nu


@nb.jit(
    nb.types.containers.Tuple((
        nb.float64,
        nb.types.EnumMember(DivType, nb.int64)
    ))(nb.complex128),
    nopython=True
)
def divergence_count(z):
    """Return a smoothed divergence count and the type of divergence."""

    delta_im = -1
    delta_re = -1
    cycle = 0
    r0 = -1
    r_drops = 0  # Counter for number of consecutive times abs(r) decreases
    a, b = z.real, z.imag
    a_cycle, b_cycle = a, b
    cutoff_re_squared = CUTOFF_RE*CUTOFF_RE

    for k in range(MAX_ITER):

        e = 0.5*np.exp(-np.pi*b)

        cycle += 1
        if cycle == 8:
            cycle = 0
            r = e*np.hypot(0.5 + a, b)/(1e-6 + np.abs(b))

            if r < r0 < 0.5:
                r_drops += 1
            else:
                r_drops = 0
            # Stop early due to likely slow linear divergence
            if r_drops >= MIN_R_DROPS:
                delta = 0.25*(CUTOFF_RE - a)
                return k + delta, DivType.LINEAR

            # Detected slow growth (maximum iterations will be reached)
            if ((k >= MIN_ITER_SLOW) and (r0 <= r)
                    and (r + (r - r0)*(MAX_ITER - k) < 8*0.05)):
                delta = 0.25*(CUTOFF_RE - a)
                return k + delta, DivType.SLOW
            r0 = r

            # Cycled back to the same value after 8 iterations
            if (a - a_cycle)**2 + (b - b_cycle)**2 < 1e-16:
                return k, DivType.CYCLE
            a_cycle = a
            b_cycle = b

        a0 = a
        # b0 = b
        s = np.sin(np.pi*a)
        c = np.cos(np.pi*a)
        # Equivalent to:
        # z = 0.25 + z - (0.25 + 0.5*z)*np.exp(np.pi*z*1j)
        # where z = a + b*1j
        a += e*(b*s - (0.5 + a)*c) + 0.25
        b -= e*(b*c + (0.5 + a0)*s)

        if b < -CUTOFF_IM:
            delta_im = get_delta_im(-b)
        if a*a + b*b > cutoff_re_squared:
            delta_re = get_delta_re(np.hypot(a, b), e)
        # Diverged by exceeding a cutoff
        if delta_im >= 0 or delta_re >= 0:
            if delta_re < 0 or delta_im <= delta_re:
                return k + delta_im, DivType.CUTOFF_IM
            else:
                return k + delta_re, DivType.CUTOFF_RE

    # Maximum iterations reached
    return -1, DivType.MAX_ITER


@nb.jit(nb.complex128(nb.float64, nb.float64), nopython=True)
def pixel_to_z(a, b):
    re = X_MIN + (X_MAX - X_MIN)*a/WIDTH
    im = Y_MAX - (Y_MAX - Y_MIN)*b/HEIGHT
    return re + 1j*im


@nb.jit(nb.float64(nb.float64), nopython=True)
def cyclic_map(g):
    """A continuous function that cycles back and forth between 0 and 1."""

    # This can be any continuous function.
    # Log scale removes high-frequency color cycles.
    # freq_div = 1
    # g = np.log1p(np.fmax(0, g/freq_div)) - np.log1p(1/freq_div)

    # Beyond this value for float64, decimals are truncated
    if g >= 2**51:
        return -1

    # Normalize and cycle
    # g += 0.5  # phase from 0 to 1
    return 1 - np.abs(2*(g - np.floor(g)) - 1)


def get_pixel(px, py):
    z = pixel_to_z(px, py)
    dc, div_type = divergence_count(z)
    match div_type:
        case DivType.MAX_ITER

授權條款

我,本作品的著作權持有者,決定用以下授權條款發佈本作品:
Creative Commons CC-Zero 此檔案在創用CC CC0 1.0 通用公有領域貢獻宣告之下分發。
在此宣告之下分發本作品者,已依據各國著作權法,在全世界放棄其對本作品所擁有的著作權及所有相關相似的法律權利,從而將本作品貢獻至公有領域。您可以複製、修改、分發和演示該作品,用於任何商業用途,所有這些都不需要請求授權。

說明

添加單行說明來描述出檔案所代表的內容
An exponential Collatz fractal with smooth coloring based on divergence speed.

在此檔案描寫的項目

描繪內容

分形 中文 (已轉換拼寫)

沒有維基數據項目的某些值

作者姓名字串 繁體中文 (已轉換拼寫):​Hugo Spinelli
維基媒體使用者名稱 繁體中文 (已轉換拼寫):​Hugo Spinelli

著作權狀態 繁體中文 (已轉換拼寫)

檔案來源 Chinese (Taiwan) (已轉換拼寫)

上傳者的原創作品 繁體中文 (已轉換拼寫)

檔案歷史

點選日期/時間以檢視該時間的檔案版本。

日期/時間縮⁠圖尺寸用戶備⁠註
目前2023年10月6日 (五) 19:37於 2023年10月6日 (五) 19:37 版本的縮圖30,720 × 17,280(41.13 MB)Hugo SpinelliUploaded own work with UploadWizard

下列頁面有用到此檔案:

全域檔案使用狀況

以下其他 wiki 使用了這個檔案: