Source code for fractalshades.models.burning_ship

# -*- coding: utf-8 -*-
import typing
import enum

import numpy as np
import mpmath
import numba

import fractalshades as fs
import fractalshades.mpmath_utils.FP_loop as fsFP

@numba.njit
def sgn(x):
    # Sign with no 0 case (0 mapped to 1)
    if x < 0.:
        return -1.
    return 1.

@numba.njit
def diffabs(X, x):
    # Evaluates |X + x| - |X|
    if X >= 0.:
        if (X + x) >= 0.:
            # /!\ Record types need to be unified: hence the * 1.
            return 1. * x
        else:
            return -(2. * X + x)
    else:
        if (X + x) <= 0.:
            return -x
        else:
            return (2. * X + x)

@numba.njit
def ddiffabsdX(X, x):
    # Evaluates |X + x| - |X|
    if X >= 0.:
        if (X + x) >= 0.:
            return 0.
        else:
            return -2.
    else:
        if (X + x) <= 0.:
            return 0.
        else:
            return 2.

@numba.njit
def ddiffabsdx(X, x):
    # Evaluates |X + x| - |X|
    if X >= 0.:
        if (X + x) >= 0.:
            return 1.
        else:
            return -1.
    else:
        if (X + x) <= 0.:
            return -1.
        else:
            return 1.


BS_flavor_list = (
    "Burning ship",
    "Perpendicular burning ship",
    "Shark fin",
    "Celtic",
    "Buffalo"
)

BS_flavor_enum =  enum.Enum(
    "BS_flavor_enum",
    BS_flavor_list,
    module=__name__
)

def get_flavor_int(flavor: str):
    # Using auto with IntEnum results in integers of increasing value,
    # starting with 1.
    return getattr(BS_flavor_enum, flavor).value

def BS_iterate(flavor_int):
    if flavor_int == 1:
        @numba.njit
        def numba_impl(xn, yn, a, b):
            # Burning ship
            return (
                xn ** 2 - yn ** 2 + a,
                2. * np.abs(xn * yn) - b
            )
    elif flavor_int == 2:
        @numba.njit
        def numba_impl(xn, yn, a, b):
            # Perpendicular burning ship
            return (
                xn ** 2 - yn ** 2 + a,
                2. * xn * np.abs(yn) - b
            )
    elif flavor_int == 3:
        @numba.njit
        def numba_impl(xn, yn, a, b):
            # Shark Fin
            return (
                xn ** 2 - yn * np.abs(yn) + a,
                2. * xn * yn - b
            )
    elif flavor_int == 4:
        @numba.njit
        def numba_impl(xn, yn, a, b):
            # Celtic
            return (
                np.abs(xn ** 2 - yn ** 2) + a,
                2. * xn * yn - b
            )
    elif flavor_int == 5:
        @numba.njit
        def numba_impl(xn, yn, a, b):
            # Buffalo
            return (
                np.abs(xn ** 2 - yn ** 2) + a,
                2. * np.abs(xn * yn) - b
            )
    return numba_impl


[docs] class Burning_ship(fs.Fractal):
[docs] def __init__( self, directory: str, flavor: typing.Literal[BS_flavor_enum]= "Burning ship" ): """ A Burning Ship Fractal (power 2). The basic equation a detailed below: .. math:: x_0 &= 0 \\\\ y_0 &= 0 \\\\ x_{n+1} &= x_n^2 - y_n^2 + a \\\\ y_{n+1} &= 2 |x_n y_n| - b where: .. math:: z_n &= x_n + i y_n \\\\ c &= a + i b Parameters ========== directory : str Path for the working base directory flavor : str The variant of Burning Ship detailed implementation, defaults to "Burning Ship". Acceptable values are listed below in notes. Notes ===== .. note:: Several variants (`flavor` parameter) are implemented with the following iteration formula: - "Perpendicular burning ship" variant of the Burning Ship Fractal. .. math:: x_{n+1} &= x_n^2 - y_n^2 + a \\\\ y_{n+1} &= 2 x_n |y_n| - b - "Shark fin" variant .. math:: x_{n+1} &= x_n^2 - y_n |y_n| + a \\\\ y_{n+1} &= 2 x_n y_n - b - "Celtic" variant .. math:: x_{n+1} &= |x_n^2 - y_n^2| + a \\\\ y_{n+1} &= 2 x_n y_n - b - "Buffalo" variant .. math:: x_{n+1} &= |x_n^2 - y_n^2| + a \\\\ y_{n+1} &= 2 |x_n y_n| - b """ super().__init__(directory) self.flavor = flavor flavor_int = get_flavor_int(flavor) # default values used for postprocessing (potential) self.potential_kind = "infinity" self.potential_d = 2 self.potential_a_d = 1. # Minimum M for valid potential estimation dic_cutoff = { "Burning ship": 1000., "Perpendicular burning ship": 1000, "Shark fin": 1.e20, "Celtic": 1000., "Buffalo": 1000, } self.potential_M_cutoff = dic_cutoff[flavor] # GUI 'badges' self.holomorphic = False self.implements_dzndc = "always" self.implements_fieldlines = True self.implements_newton = False self.implements_Milnor = False self.implements_interior_detection = "no" self.implements_deepzoom = False self.xnyn_iterate = BS_iterate(flavor_int) # Cache for numba compiled implementations self._numba_cache = {} self._numba_iterate_cache = { "calc_std_div": ((), None), } # args, numba_impl
[docs] @fs.utils.calc_options def calc_std_div(self, *, calc_name: str, subset, max_iter: int, M_divergence: float, calc_orbit: bool = False, backshift: int = 0 ): """ Basic iterations for Burning ship set. Parameters ========== calc_name : str The string identifier for this calculation subset : `fractalshades.postproc.Fractal_array` A boolean array-like, where False no calculation is performed If `None`, all points are calculated. Defaults to `None`. max_iter : int the maximum iteration number. If reached, the loop is exited with exit code "max_iter". M_divergence : float The diverging radius. If reached, the loop is exited with exit code "divergence" calc_orbit: bool If True, stores the value of an orbit point @ exit - orbit_shift backshift: int (> 0) The count of iterations backward for the stored orbit starting point Notes ===== The following complex fields will be calculated: *xn* *yn* and its derivatives (*dxnda*, *dxndb*, *dynda*, *dyndb*). If calc_orbit is activated, *zn_orbit* will also be stored Exit codes are *max_iter*, *divergence*. """ complex_codes = ["xn", "yn", "dxnda", "dxndb", "dynda", "dyndb"] xn = 0 yn = 1 dxnda = 2 dxndb = 3 dynda = 4 dyndb = 5 i_xnorbit = -1 i_ynorbit = -1 if calc_orbit: complex_codes += ["xn_orbit", "yn_orbit"] i_xnorbit = 6 i_ynorbit = 7 int_codes = [] stop_codes = ["max_iter", "divergence", "stationnary"] def set_state(): def impl(instance): instance.codes = (complex_codes, int_codes, stop_codes) instance.complex_type = np.float64 instance.potential_M = M_divergence if calc_orbit: instance.backshift = backshift else: instance.backshift = None return impl def initialize(): @numba.njit def numba_init_impl(Z, U, c): # Not much to do here pass return numba_init_impl Mdiv_sq = M_divergence ** 2 flavor_int = getattr(BS_flavor_enum, self.flavor).value iterate_once = self.from_numba_cache( "iterate_once", flavor_int, xn, yn, dxnda, dxndb, dynda, dyndb, max_iter, Mdiv_sq ) xnyn_iterate = self.xnyn_iterate def iterate(): new_args = ( calc_orbit, i_xnorbit, i_ynorbit, backshift, xn, yn, iterate_once, xnyn_iterate ) args, numba_impl = self._numba_iterate_cache["calc_std_div"] if new_args == args: return numba_impl numba_impl = fs.core.numba_iterate_BS(*new_args) self._numba_iterate_cache["calc_std_div"] = (new_args, numba_impl) return numba_impl return { "set_state": set_state, "initialize": initialize, "iterate": iterate }
def from_numba_cache( self, key, flavor_int, xn, yn, dxnda, dxndb, dynda, dyndb, max_iter, Mdiv_sq ): """ Returns the numba implementation if exists to avoid unnecessary recompilation""" cache = self._numba_cache full_key = (key, flavor_int, xn, yn, dxnda, dxndb, dynda, dyndb, max_iter, Mdiv_sq) try: return cache[full_key] except KeyError: if key == "iterate_once": @numba.njit def numba_impl(c, Z, U, stop_reason, n_iter): if n_iter >= max_iter: stop_reason[0] = 0 return 1 a = c.real b = c.imag X = Z[xn] Y = Z[yn] dXdA = Z[dxnda] dXdB = Z[dxndb] dYdA = Z[dynda] dYdB = Z[dyndb] if flavor_int == 1: # Standard BS implementation Z[xn] = X ** 2 - Y ** 2 + a Z[yn] = 2. * np.abs(X * Y) - b # Jacobian Z[dxnda] = 2 * (X * dXdA - Y * dYdA) + 1. Z[dxndb] = 2 * (X * dXdB - Y * dYdB) Z[dynda] = 2 * (np.abs(X) * sgn(Y) * dYdA + sgn(X) * dXdA * np.abs(Y)) Z[dyndb] = 2 * (np.abs(X) * sgn(Y) * dYdB + sgn(X) * dXdB * np.abs(Y)) - 1. elif flavor_int == 2: # Perpendicular BS implementation Z[xn] = X ** 2 - Y ** 2 + a Z[yn] = 2. * X * np.abs(Y) - b # Jacobian Z[dxnda] = 2 * (X * dXdA - Y * dYdA) + 1. Z[dxndb] = 2 * (X * dXdB - Y * dYdB) Z[dynda] = 2 * (X * sgn(Y) * dYdA + dXdA * np.abs(Y)) Z[dyndb] = 2 * (X * sgn(Y) * dYdB + dXdB * np.abs(Y)) - 1. elif flavor_int == 3: # Shark Fin Z[xn] = X ** 2 - Y * np.abs(Y) + a Z[yn] = 2. * X * Y - b # Jacobian Z[dxnda] = 2 * (X * dXdA - np.abs(Y) * dYdA) + 1. Z[dxndb] = 2 * (X * dXdB - np.abs(Y) * dYdB) Z[dynda] = 2. * (dXdA * Y + X * dYdA) Z[dyndb] = 2. * (dXdB * Y + X * dYdB) - 1. elif flavor_int == 4: # Celtic x2my2 = X ** 2 - Y ** 2 Z[xn] = np.abs(x2my2) + a Z[yn] = 2. * X * Y - b # Jacobian Z[dxnda] = 2. * sgn(x2my2) * (X * dXdA - Y * dYdA) Z[dxndb] = 2. * sgn(x2my2) * (X * dXdB - Y * dYdB) Z[dynda] = 2. * (dXdA * Y + X * dYdA) Z[dyndb] = 2. * (dXdB * Y + X * dYdB) - 1. elif flavor_int == 5: # Buffalo x2my2 = X ** 2 - Y ** 2 Z[xn] = np.abs(x2my2) + a Z[yn] = 2. * np.abs(X * Y) - b # Jacobian Z[dxnda] = 2. * sgn(x2my2) * (X * dXdA - Y * dYdA) Z[dxndb] = 2. * sgn(x2my2) * (X * dXdB - Y * dYdB) Z[dynda] = 2. * (np.abs(X) * sgn(Y) * dYdA + sgn(X) * dXdA * np.abs(Y)) Z[dyndb] = 2. * (np.abs(X) * sgn(Y) * dYdB + sgn(X) * dXdB * np.abs(Y)) - 1. if Z[xn] ** 2 + Z[yn] ** 2 > Mdiv_sq: stop_reason[0] = 1 return 1 return 0 else: raise NotImplementedError(key) cache[full_key] = numba_impl return numba_impl @fs.utils.interactive_options def coords(self, x, y, pix, dps): return super().coords(x, y, pix, dps)
#============================================================================== #============================================================================== # Utility functions used in perturbation burning ship def _dfxdx(flavor_int): if flavor_int == 1: @numba.njit def numba_impl(x, y): return 2. * x elif flavor_int == 2: @numba.njit def numba_impl(x, y): return 2. * x elif flavor_int == 3: @numba.njit def numba_impl(x, y): return 2. * x elif flavor_int == 4 or flavor_int == 5: # x_{n+1} &= |x_n^2 - y_n^2| + a \\\\ @numba.njit def numba_impl(x, y): x2my2 = x * x - y * y return 2. * sgn(x2my2) * x return numba_impl def _dfxdy(flavor_int): if flavor_int == 1: @numba.njit def numba_impl(x, y): return -2. * y elif flavor_int == 2: @numba.njit def numba_impl(x, y): return -2. * y elif flavor_int == 3: @numba.njit def numba_impl(x, y): return -2. * np.abs(y) elif flavor_int == 4 or flavor_int == 5: # x_{n+1} &= |x_n^2 - y_n^2| + a \\\\ @numba.njit def numba_impl(x, y): x2my2 = x * x - y * y return -2. * sgn(x2my2) * y return numba_impl def _dfydx(flavor_int): if flavor_int == 1: @numba.njit def numba_impl(x, y): return 2. * sgn(x) * np.abs(y) elif flavor_int == 2: @numba.njit def numba_impl(x, y): return 2. * np.abs(y) elif flavor_int == 3: @numba.njit def numba_impl(x, y): return 2. * y elif flavor_int == 4: @numba.njit def numba_impl(x, y): return 2. * y elif flavor_int == 5: @numba.njit def numba_impl(x, y): return 2. * sgn(x) * np.abs(y) return numba_impl def _dfydy(flavor_int): if flavor_int == 1: @numba.njit def numba_impl(x, y): return 2. * sgn(y) * np.abs(x) elif flavor_int == 2: @numba.njit def numba_impl(x, y): return 2. * sgn(y) * x elif flavor_int == 3: @numba.njit def numba_impl(x, y): return 2. * x elif flavor_int == 4: @numba.njit def numba_impl(x, y): return 2. * x elif flavor_int == 5: @numba.njit def numba_impl(x, y): return 2. * sgn(y) * np.abs(x) return numba_impl def _p_iter_zn(flavor_int, xn, yn): """ Perturbation iteration for different burning ship 'flavors'""" if flavor_int == 1: @numba.njit # Burning Ship def numba_impl(Z, ref_xn, ref_yn, a, b): # Modifies in-place xn, yn ref_xyn = ref_xn * ref_yn new_xn = ( Z[xn] * (Z[xn] + 2. * ref_xn) - Z[yn] * (Z[yn] + 2. * ref_yn) + a ) new_yn = ( 2. * diffabs( ref_xyn, Z[xn] * Z[yn] + Z[xn] * ref_yn + Z[yn] * ref_xn ) - b ) Z[xn] = new_xn Z[yn] = new_yn elif flavor_int == 2: # Perpendicular BS @numba.njit def numba_impl(Z, ref_xn, ref_yn, a, b): new_xn = ( Z[xn] * (Z[xn] + 2. * ref_xn) - Z[yn] * (Z[yn] + 2. * ref_yn) + a ) new_yn = ( 2. * ( ref_xn * diffabs(ref_yn, Z[yn]) + Z[xn] * np.abs(ref_yn + Z[yn]) ) - b ) Z[xn] = new_xn Z[yn] = new_yn elif flavor_int == 3: # Shark Fin @numba.njit def numba_impl(Z, ref_xn, ref_yn, a, b): new_xn = ( Z[xn] * (Z[xn] + 2. * ref_xn) - ref_yn * diffabs(ref_yn, Z[yn]) - Z[yn] * np.abs(ref_yn + Z[yn]) + a ) new_yn = 2. * (ref_xn * Z[yn] + ref_yn * Z[xn] + Z[xn] * Z[yn]) - b Z[xn] = new_xn Z[yn] = new_yn elif flavor_int == 4: # Celtic @numba.njit def numba_impl(Z, ref_xn, ref_yn, a, b): ref_x2my2 = ref_xn * ref_xn - ref_yn * ref_yn new_xn = ( diffabs( ref_x2my2, Z[xn] * (Z[xn] + 2. * ref_xn) - Z[yn] * (Z[yn] + 2. * ref_yn) ) + a ) new_yn = 2. * (ref_xn * Z[yn] + ref_yn * Z[xn] + Z[xn] * Z[yn]) - b Z[xn] = new_xn Z[yn] = new_yn elif flavor_int == 5: # Buffalo @numba.njit def numba_impl(Z, ref_xn, ref_yn, a, b): ref_xyn = ref_xn * ref_yn ref_x2my2 = ref_xn * ref_xn - ref_yn * ref_yn new_xn = ( diffabs( ref_x2my2, Z[xn] * (Z[xn] + 2. * ref_xn) - Z[yn] * (Z[yn] + 2. * ref_yn) ) + a ) new_yn = ( 2. * diffabs( ref_xyn, Z[xn] * Z[yn] + Z[xn] * ref_yn + Z[yn] * ref_xn ) - b ) Z[xn] = new_xn Z[yn] = new_yn return numba_impl def _p_iter_hessian(flavor_int, xn, yn, dxnda, dxndb, dynda, dyndb): """ Hessian matrix perturbation iteration for different Burning ship 'flavors' https://fractalforums.org/fractal-mathematics-and-new-theories/28/perturbation-theory/487/msg3226#msg3226 """ if flavor_int == 1: # Burning Ship @numba.njit def numba_impl( Z, ref_xn, ref_yn, ref_dxnda, ref_dxndb, ref_dynda, ref_dyndb ): # Modifies in-place the Hessian matrix ref_xyn = ref_xn * ref_yn _opX = ref_xyn d_opX_da = ref_dxnda * ref_yn + ref_xn * ref_dynda d_opX_db = ref_dxndb * ref_yn + ref_xn * ref_dyndb _opx = Z[xn] * Z[yn] + Z[xn] * ref_yn + Z[yn] * ref_xn d_opx_da = ( Z[dxnda] * Z[yn] + Z[xn] * Z[dynda] + Z[dxnda] * ref_yn + Z[xn] * ref_dynda + Z[dynda] * ref_xn + Z[yn] * ref_dxnda ) d_opx_db = ( Z[dxndb] * Z[yn] + Z[xn] * Z[dyndb] + Z[dxndb] * ref_yn + Z[xn] * ref_dyndb + Z[dyndb] * ref_xn + Z[yn] * ref_dxndb ) _ddiffabsdX = ddiffabsdX(_opX, _opx) _ddiffabsdx = ddiffabsdx(_opX, _opx) new_dxnda = ( 2. * ((ref_xn + Z[xn]) * Z[dxnda] + ref_dxnda * Z[xn]) -2. * ((ref_yn + Z[yn]) * Z[dynda] + ref_dynda * Z[yn]) ) new_dxndb = ( 2. * ((ref_xn + Z[xn]) * Z[dxndb] + ref_dxndb * Z[xn]) -2. * ((ref_yn + Z[yn]) * Z[dyndb] + ref_dyndb * Z[yn]) ) new_dynda = 2. * (_ddiffabsdX * d_opX_da + _ddiffabsdx * d_opx_da) new_dyndb = 2. * (_ddiffabsdX * d_opX_db + _ddiffabsdx * d_opx_db) Z[dxnda] = new_dxnda Z[dxndb] = new_dxndb Z[dynda] = new_dynda Z[dyndb] = new_dyndb elif flavor_int == 2: # Perpendicular BS @numba.njit def numba_impl( Z, ref_xn, ref_yn, ref_dxnda, ref_dxndb, ref_dynda, ref_dyndb ): _diffabs = diffabs(ref_yn, Z[yn]) _ddiffabsdX = ddiffabsdX(ref_yn, Z[yn]) _ddiffabsdx = ddiffabsdx(ref_yn, Z[yn]) Y_y = ref_yn + Z[yn] _abs = np.abs(Y_y) _sgn = sgn(Y_y) new_dxnda = 2. * ( ((ref_xn + Z[xn]) * Z[dxnda] + ref_dxnda * Z[xn]) -((ref_yn + Z[yn]) * Z[dynda] + ref_dynda * Z[yn]) ) new_dxndb = 2. * ( ((ref_xn + Z[xn]) * Z[dxndb] + ref_dxndb * Z[xn]) -((ref_yn + Z[yn]) * Z[dyndb] + ref_dyndb * Z[yn]) ) new_dynda = 2. * ( ref_dxnda * _diffabs + ref_xn * ( _ddiffabsdX * ref_dynda + _ddiffabsdx * Z[dynda] ) + Z[dxnda] * _abs + Z[xn] * _sgn * (ref_dynda + Z[dynda]) ) new_dyndb = 2. * ( ref_dxndb * _diffabs + ref_xn * ( _ddiffabsdX * ref_dyndb + _ddiffabsdx * Z[dyndb] ) + Z[dxndb] * _abs + Z[xn] * _sgn * (ref_dyndb + Z[dyndb]) ) Z[dxnda] = new_dxnda Z[dxndb] = new_dxndb Z[dynda] = new_dynda Z[dyndb] = new_dyndb elif flavor_int == 3: # Shark Fin @numba.njit def numba_impl( Z, ref_xn, ref_yn, ref_dxnda, ref_dxndb, ref_dynda, ref_dyndb ): _diffabs = diffabs(ref_yn, Z[yn]) _ddiffabsdX = ddiffabsdX(ref_yn, Z[yn]) _ddiffabsdx = ddiffabsdx(ref_yn, Z[yn]) Y_y = ref_yn + Z[yn] _abs = np.abs(Y_y) _sgn = sgn(Y_y) new_dxnda = ( Z[dxnda] * (Z[xn] + 2. * ref_xn) + Z[xn] * (Z[dxnda] + 2. * ref_dxnda) - ref_dynda * _diffabs - ref_yn * (ref_dynda * _ddiffabsdX + Z[dynda] * _ddiffabsdx) - Z[dynda] * _abs - Z[yn] * _sgn * (ref_dynda + Z[dynda]) ) new_dxndb = ( Z[dxndb] * (Z[xn] + 2. * ref_xn) + Z[xn] * (Z[dxndb] + 2. * ref_dxndb) - ref_dyndb * _diffabs - ref_yn * (ref_dyndb * _ddiffabsdX + Z[dyndb] * _ddiffabsdx) - Z[dyndb] * _abs - Z[yn] * _sgn * (ref_dyndb + Z[dyndb]) ) new_dynda = 2. * ( ref_dxnda * Z[yn] + ref_xn * Z[dynda] + ref_dynda * Z[xn] + ref_yn * Z[dxnda] + Z[dxnda] * Z[yn] + Z[xn] * Z[dynda] ) new_dyndb = 2. * ( ref_dxndb * Z[yn] + ref_xn * Z[dyndb] + ref_dyndb * Z[xn] + ref_yn * Z[dxndb] + Z[dxndb] * Z[yn] + Z[xn] * Z[dyndb] ) Z[dxnda] = new_dxnda Z[dxndb] = new_dxndb Z[dynda] = new_dynda Z[dyndb] = new_dyndb elif flavor_int == 4: # Celtic @numba.njit def numba_impl( Z, ref_xn, ref_yn, ref_dxnda, ref_dxndb, ref_dynda, ref_dyndb ): # Modifies in-place the Hessian matrix # Implementation for dxnda, dxndb # new_xn = ( # diffabs( # ref_x2my2, # Z[xn] * (Z[xn] + 2. * ref_xn) - Z[yn] * (Z[yn] + 2. * ref_yn) # ) + a # ) ref_x2my2 = ref_xn * ref_xn - ref_yn * ref_yn _opX = ref_x2my2 d_opX_da = 2. * (ref_xn * ref_dxnda - ref_yn * ref_dynda) d_opX_db = 2. * (ref_xn * ref_dxndb - ref_yn * ref_dyndb) _opx = Z[xn] * (Z[xn] + 2. * ref_xn) - Z[yn] * (Z[yn] + 2. * ref_yn) d_opx_da = ( Z[dxnda] * (Z[xn] + 2. * ref_xn) + Z[xn] * (Z[dxnda] + 2. * ref_dxnda) - Z[dynda] * (Z[yn] + 2. * ref_yn) - Z[yn] * (Z[dynda] + 2. * ref_dynda) ) d_opx_db = ( Z[dxndb] * (Z[xn] + 2. * ref_xn) + Z[xn] * (Z[dxndb] + 2. * ref_dxndb) - Z[dyndb] * (Z[yn] + 2. * ref_yn) - Z[yn] * (Z[dyndb] + 2. * ref_dyndb) ) _ddiffabsdX = ddiffabsdX(_opX, _opx) _ddiffabsdx = ddiffabsdx(_opX, _opx) new_dxnda = _ddiffabsdX * d_opX_da + _ddiffabsdx * d_opx_da new_dxndb = _ddiffabsdX * d_opX_db + _ddiffabsdx * d_opx_db # Implementation based on Shark Fin for dynda, dyndb new_dynda = 2. * ( ref_dxnda * Z[yn] + ref_xn * Z[dynda] + ref_dynda * Z[xn] + ref_yn * Z[dxnda] + Z[dxnda] * Z[yn] + Z[xn] * Z[dynda] ) new_dyndb = 2. * ( ref_dxndb * Z[yn] + ref_xn * Z[dyndb] + ref_dyndb * Z[xn] + ref_yn * Z[dxndb] + Z[dxndb] * Z[yn] + Z[xn] * Z[dyndb] ) Z[dxnda] = new_dxnda Z[dxndb] = new_dxndb Z[dynda] = new_dynda Z[dyndb] = new_dyndb elif flavor_int == 5: # Buffalo @numba.njit def numba_impl( Z, ref_xn, ref_yn, ref_dxnda, ref_dxndb, ref_dynda, ref_dyndb ): # Modifies in-place the Hessian matrix # Implementation based on Celtic for dxnda, dxndb ref_x2my2 = ref_xn * ref_xn - ref_yn * ref_yn _opX = ref_x2my2 d_opX_da = 2. * (ref_xn * ref_dxnda - ref_yn * ref_dynda) d_opX_db = 2. * (ref_xn * ref_dxndb - ref_yn * ref_dyndb) _opx = Z[xn] * (Z[xn] + 2. * ref_xn) - Z[yn] * (Z[yn] + 2. * ref_yn) d_opx_da = ( Z[dxnda] * (Z[xn] + 2. * ref_xn) + Z[xn] * (Z[dxnda] + 2. * ref_dxnda) - Z[dynda] * (Z[yn] + 2. * ref_yn) - Z[yn] * (Z[dynda] + 2. * ref_dynda) ) d_opx_db = ( Z[dxndb] * (Z[xn] + 2. * ref_xn) + Z[xn] * (Z[dxndb] + 2. * ref_dxndb) - Z[dyndb] * (Z[yn] + 2. * ref_yn) - Z[yn] * (Z[dyndb] + 2. * ref_dyndb) ) _ddiffabsdX = ddiffabsdX(_opX, _opx) _ddiffabsdx = ddiffabsdx(_opX, _opx) new_dxnda = _ddiffabsdX * d_opX_da + _ddiffabsdx * d_opx_da new_dxndb = _ddiffabsdX * d_opX_db + _ddiffabsdx * d_opx_db # Implementation based on Burning Ship for dynda, dyndb ref_xyn = ref_xn * ref_yn _opX = ref_xyn d_opX_da = ref_dxnda * ref_yn + ref_xn * ref_dynda d_opX_db = ref_dxndb * ref_yn + ref_xn * ref_dyndb _opx = Z[xn] * Z[yn] + Z[xn] * ref_yn + Z[yn] * ref_xn d_opx_da = ( Z[dxnda] * Z[yn] + Z[xn] * Z[dynda] + Z[dxnda] * ref_yn + Z[xn] * ref_dynda + Z[dynda] * ref_xn + Z[yn] * ref_dxnda ) d_opx_db = ( Z[dxndb] * Z[yn] + Z[xn] * Z[dyndb] + Z[dxndb] * ref_yn + Z[xn] * ref_dyndb + Z[dyndb] * ref_xn + Z[yn] * ref_dxndb ) _ddiffabsdX = ddiffabsdX(_opX, _opx) _ddiffabsdx = ddiffabsdx(_opX, _opx) new_dynda = 2. * (_ddiffabsdX * d_opX_da + _ddiffabsdx * d_opx_da) new_dyndb = 2. * (_ddiffabsdX * d_opX_db + _ddiffabsdx * d_opx_db) Z[dxnda] = new_dxnda Z[dxndb] = new_dxndb Z[dynda] = new_dynda Z[dyndb] = new_dyndb return numba_impl #============================================================================== #==============================================================================
[docs] class Perturbation_burning_ship(fs.PerturbationFractal):
[docs] def __init__( self, directory: str, flavor: typing.Literal[BS_flavor_enum]= "Burning ship" ): """ Arbitrary-precision class for the Burning ship and other "abs variations". The Burning Ship fractal, first described by Michael Michelitsch and Otto E. Rössler in 1992, is a variant of the mandelbrot fractal which involve the absolute value function, making the formula non-analytic: .. math:: x_0 &= 0 \\\\ y_0 &= 0 \\\\ x_{n+1} &= x_n^2 - y_n^2 + a \\\\ y_{n+1} &= 2 |x_n y_n| - b where: .. math:: z_n &= x_n + i y_n \\\\ c &= a + i b For a more comprehensive introduction, we recommend the paper *At the Helm of the Burning Ship* :cite:p:`burning_ship`. This class implements arbitrary precision for the reference orbit, ball method period search, newton search, perturbation method, chained billinear approximations. Parameters ---------- directory: str Path for the working base directory flavor: str The variant of Burning Ship detailed implementation, defaults to "Burning Ship". Acceptable values and iteration formula are detailed in the notes of `fractalshades.models.Burning_ship` documentation. """ super().__init__(directory) self.flavor = flavor flavor_int = get_flavor_int(flavor) # Sets default values used for postprocessing (potential) self.potential_kind = "infinity" self.potential_d = 2 self.potential_a_d = 1. self.potential_M_cutoff = 1000. # Minimum M for valid potential # Set parameters for the full precision orbit self.critical_pt = 0. self.FP_code = ["xn", "yn"] self.holomorphic = False # GUI 'badges' self.holomorphic = False self.implements_dzndc = "always" self.implements_fieldlines = True self.implements_newton = False self.implements_Milnor = False self.implements_interior_detection = "no" self.implements_deepzoom = True # Orbit 'on the fly' calculation self.xnyn_iterate = BS_iterate(flavor_int) # Cache for numba compiled implementations self._numba_cache = {} self._numba_initialize_cache = ((), None) # args, numba_impl self._numba_iterate_cache = ((), None) # args, numba_impl
def FP_loop(self, NP_orbit, c0): """ The full precision loop ; fills in place NP_orbit """ xr_detect_activated = self.xr_detect_activated max_orbit_iter = (NP_orbit.shape)[0] - 1 M_divergence = self.M_divergence flavor_int = get_flavor_int(self.flavor) x = c0.real y = c0.imag seed_prec = mpmath.mp.prec (i, partial_dict, xr_dict ) = fsFP.perturbation_nonholomorphic_FP_loop( NP_orbit.view(dtype=np.float64), xr_detect_activated, max_orbit_iter, M_divergence * 2, # to be sure ref exit after close points str(x).encode('utf8'), str(y).encode('utf8'), seed_prec, kind=flavor_int ) return i, partial_dict, xr_dict
[docs] @fs.utils.calc_options def calc_std_div(self, *, calc_name: str, subset, max_iter: int, M_divergence: float, BLA_eps: float = 1e-6, calc_hessian: bool = True, calc_orbit: bool = False, backshift: int = 0 ): """ Perturbation iterations (arbitrary precision) for this class. Parameters ========== calc_name : str The string identifier for this calculation subset : A boolean array-like, where False no calculation is performed If `None`, all points are calculated. Defaults to `None`. max_iter : int the maximum iteration number. If reached, the loop is exited with exit code "max_iter". M_divergence : float The diverging radius. If reached, the loop is exited with exit code "divergence" BLA_eps: float Relative error criteria for BLA (default: 1.e-6) If `None`, BLA is not activated. calc_hessian: bool if True, the derivatives will be caculated allowing distance estimation and shading. calc_orbit: bool If True, stores the value of an orbit point @ exit - orbit_shift backshift: int (> 0) The count of iterations backward for the stored orbit starting point """ complex_codes = ["xn", "yn"] xn = 0 yn = 1 code_int = 2 if calc_hessian: complex_codes += ["dxnda", "dxndb", "dynda", "dyndb"] dxnda = code_int + 0 dxndb = code_int + 1 dynda = code_int + 2 dyndb = code_int + 3 code_int += 4 else: dxnda = -1 dxndb = -1 dynda = -1 dyndb = -1 if calc_orbit: complex_codes += ["xn_orbit", "yn_orbit"] i_xnorbit = code_int + 0 i_ynorbit = code_int + 1 code_int += 2 else: i_xnorbit = i_ynorbit = -1 # Integer int32 fields codes "U" int_codes = ["ref_cycle_iter"] # Position in ref orbit # Stop codes stop_codes = ["max_iter", "divergence"] reason_max_iter = 0 reason_M_divergence = 1 #---------------------------------------------------------------------- # Define the functions used for BLA approximation # BLA triggered ? BLA_activated = ( (BLA_eps is not None) and (self.dx < fs.settings.newton_zoom_level) ) # H = [dfxdx dfxdy] [dfx] = H x [dx] # [dfydx dfydy] [dfy] [dy] flavor_int = get_flavor_int(self.flavor) cache_args = (flavor_int, xn, yn, dxnda, dxndb, dynda, dyndb) dfxdx = self.from_numba_cache("dfxdx", *cache_args) dfxdy = self.from_numba_cache("dfxdy", *cache_args) dfydx = self.from_numba_cache("dfydx", *cache_args) dfydy = self.from_numba_cache("dfydy", *cache_args) def set_state(): def impl(instance): instance.complex_type = np.float64 instance.potential_M = M_divergence instance.codes = (complex_codes, int_codes, stop_codes) instance.dfxdx = dfxdx # (flavor_int) instance.dfxdy = dfxdy # (flavor_int) instance.dfydx = dfydx # (flavor_int) instance.dfydy = dfydy # (flavor_int) if calc_orbit: instance.backshift = backshift else: instance.backshift = None return impl #---------------------------------------------------------------------- # Defines initialize - jitted implementation def initialize(): new_args = (xn, yn, dxnda, dxndb, dynda, dyndb) args, numba_impl = self._numba_initialize_cache if new_args == args: return numba_impl numba_impl = fs.perturbation.numba_initialize_BS(*new_args) self._numba_initialize_cache = (new_args, numba_impl) return numba_impl #---------------------------------------------------------------------- # Defines iterate - jitted implementation M_divergence_sq = M_divergence ** 2 # Xr triggered for ultra-deep zoom xr_detect_activated = self.xr_detect_activated p_iter_zn = self.from_numba_cache("p_iter_zn", *cache_args) p_iter_hessian = self.from_numba_cache("p_iter_hessian", *cache_args) xnyn_iterate = self.xnyn_iterate def iterate(): new_args = ( M_divergence_sq, max_iter, reason_max_iter, reason_M_divergence, xr_detect_activated, BLA_activated, calc_hessian, xn, yn, dxnda, dxndb, dynda, dyndb, p_iter_zn, p_iter_hessian, calc_orbit, i_xnorbit, i_ynorbit, backshift, xnyn_iterate ) args, numba_impl = self._numba_iterate_cache if new_args == args: return numba_impl numba_impl = fs.perturbation.numba_iterate_BS(*new_args) self._numba_iterate_cache = (new_args, numba_impl) return numba_impl return { "set_state": set_state, "initialize": initialize, "iterate": iterate }
def from_numba_cache(self, key, flavor_int, xn, yn, dxnda, dxndb, dynda, dyndb): """ Returns the numba implementation if exists to avoid unnecessary recompilation""" cache = self._numba_cache full_key = (key, flavor_int, xn, yn, dxnda, dxndb, dynda, dyndb) try: return cache[full_key] except KeyError: if key == "dfxdx": numba_impl = _dfxdx(flavor_int) elif key == "dfxdy": numba_impl = _dfxdy(flavor_int) elif key == "dfydx": numba_impl = _dfydx(flavor_int) elif key == "dfydy": numba_impl = _dfydy(flavor_int) elif key == "p_iter_zn": numba_impl = _p_iter_zn(flavor_int, xn, yn) elif key == "p_iter_hessian": numba_impl = _p_iter_hessian( flavor_int, xn, yn, dxnda, dxndb, dynda, dyndb ) else: raise NotImplementedError(key) cache[full_key] = numba_impl return numba_impl #------------------------------------------------------------------------------ # Newton search & other related methods def _ball_method(self, c, px, maxiter, M_divergence): """ Order 1 ball method: Cython wrapper""" x = c.real y = c.imag seed_prec = mpmath.mp.prec flavor_int = get_flavor_int(self.flavor) order = fsFP.perturbation_nonholomorphic_ball_method( str(x).encode('utf8'), str(y).encode('utf8'), seed_prec, str(px).encode('utf8'), maxiter, M_divergence, kind=flavor_int ) if order == -1: return None return order @staticmethod def find_nucleus(c, order, eps_pixel, max_newton=None, eps_cv=None): """ Run Newton search to find z0 so that f^n(z0) == 0 : Cython wrapper Includes a "divide by undesired roots" technique so that solutions using divisors of n are disregarded. # eps_pixel = self.dx * (1. / self.nx) """ raise NotImplementedError("Divide by undesired roots technique " "not implemented for burning ship") def find_any_nucleus( self, c, order, eps_pixel, max_newton=None, eps_cv=None ): """ Run Newton search to find z0 so that f^n(z0) == 0 : Cython wrapper """ if order is None: raise ValueError("order shall be defined for Newton method") x = c.real y = c.imag seed_prec = mpmath.mp.prec if max_newton is None: max_newton = 80 if eps_cv is None: eps_cv = mpmath.mpf(val=(2, -seed_prec)) flavor_int = get_flavor_int(self.flavor) is_ok, val = fsFP.perturbation_nonholomorphic_find_any_nucleus( str(x).encode('utf8'), str(y).encode('utf8'), seed_prec, order, max_newton, str(eps_cv).encode('utf8'), str(eps_pixel).encode('utf8'), kind=flavor_int ) return is_ok, val def _nucleus_size_estimate(self, c0, order): """ Nucleus size estimate Parameters ---------- c0 : position of the nucleus order : cycle order Returns ------- nucleus_size : size estimate of the nucleus julia_size : size estimate of the Julian embedded set https://mathr.co.uk/blog/2016-12-24_deriving_the_size_estimate.html Structure in the parameter dependence of order and chaos for the quadratic map Brian R Hunt and Edward Ott J. Phys. A: Math. Gen. 30 (1997) 7067–7076 https://fractalforums.org/fractal-mathematics-and-new-theories/28/miniset-and-embedded-julia-size-estimates/912/msg4805#msg4805 julia size estimate : r_J = r_M ** ((n+1)*(n-1)/n**2) """ x = c0.real y = c0.imag seed_prec = mpmath.mp.prec flavor_int = get_flavor_int(self.flavor) (nucleus_size, skew ) = fsFP.perturbation_nonholomorphic_nucleus_size_estimate( str(x).encode('utf8'), str(y).encode('utf8'), seed_prec, order, kind=flavor_int ) # r_J = r_M ** 0.75 for power 2 Mandelbrot sqrt = np.sqrt(nucleus_size) sqrtsqrt = np.sqrt(sqrt) julia_size = sqrtsqrt * sqrt return nucleus_size, julia_size, skew #============================================================================== # GUI : "interactive options" #============================================================================== @fs.utils.interactive_options def coords(self, x, y, pix, dps): return super().coords(x, y, pix, dps) @fs.utils.interactive_options def ball_method_order(self, x, y, pix, dps, maxiter: int=100000, radius_pixels: int=25): return super().ball_method_order(x, y, pix, dps, maxiter, radius_pixels) @fs.utils.interactive_options def newton_search(self, x, y, pix, dps, maxiter: int=100000, radius_pixels: int=3): ( x_str, y_str, maxiter, radius_pixels, radius_str, dps, order, xn_str, yn_str, size_estimates ) = self._newton_search( x, y, pix, dps, maxiter, radius_pixels ) if size_estimates is not None: (nucleus_size, julia_size, skew) = size_estimates else: nucleus_size = None julia_size = None skew = np.array(((np.nan, np.nan), (np.nan, np.nan))) res_str = f""" newton_search = {{ "x_start": "{x_str}", "y_start": "{y_str}", "maxiter": {maxiter}, "radius_pixels": {radius_pixels}, "radius": "{radius_str}", "calculation dps": {dps} "order": {order} "x_nucleus": "{xn_str}", "y_nucleus": "{yn_str}", "nucleus_size": "{nucleus_size}", "julia_size": "{julia_size}", "skew_00": "{skew[0, 0]}", "skew_01": "{skew[0, 1]}", "skew_10": "{skew[1, 0]}", "skew_11": "{skew[1, 1]}", }} """ return res_str @fs.utils.interactive_options def quick_skew_estimate(self, x, y, pix, dps, maxiter: int=100000): flavor_int = get_flavor_int(self.flavor) seed_prec = mpmath.mp.prec skew = fsFP.perturbation_nonholomorphic_skew_estimate( str(x).encode('utf8'), str(y).encode('utf8'), seed_prec, maxiter, 100000., kind=flavor_int ) res_str = f""" quick_skew_estimate = {{ "skew_00": "{skew[0, 0]}", "skew_01": "{skew[0, 1]}", "skew_10": "{skew[1, 0]}", "skew_11": "{skew[1, 1]}", }} """ return res_str