Source code for fractalshades.models.mandelbrot_Mn

# -*- coding: utf-8 -*-
import math
from typing import Optional

import numpy as np
import numba
import mpmath

import fractalshades as fs
import fractalshades.mpmath_utils.FP_loop as fsFP
from fractalshades.postproc import Fractal_array

def Mn_iterate(n):
    @numba.njit
    def numba_impl(zn, c):
        return zn ** n + c
    return numba_impl


[docs] class Mandelbrot_N(fs.Fractal):
[docs] def __init__(self, directory: str, exponent: int): """ A standard power-N Mandelbrot Fractal set implementation. .. math:: z_0 &= 0 \\\\ z_{n+1} &= {z_{n}}^N + c Parameters ========== directory : str Path for the working base directory """ super().__init__(directory) self.exponent = exponent # Needed for serialization # default values used for postprocessing (potential) self.potential_kind = "infinity" self.potential_d = exponent self.potential_a_d = 1. # GUI 'badges' self.holomorphic = True self.implements_dzndc = "always" self.implements_fieldlines = True self.implements_newton = True self.implements_Milnor = True self.implements_interior_detection = "always" self.implements_deepzoom = False self.zn_iterate = Mn_iterate(exponent) # Cache for numba compiled implementations self._numba_cache = {} self._numba_iterate_cache = { "calc_std_div": ((), None), "newton_calc": ((), None), } # args, numba_impl
[docs] @fs.utils.calc_options def calc_std_div(self, *, calc_name: str, subset: Optional[Fractal_array] = None, max_iter: int, M_divergence: float, epsilon_stationnary: float, calc_d2zndc2: bool = False, calc_orbit: bool = False, backshift: int = 0 ): """ Basic iterations for the Mandelbrot standard set (power n). Parameters ========== calc_name : str The string identifier for this calculation subset : Optional `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" epsilon_stationnary : float A small float to early exit non-divergent cycles (based on cumulated dzndz product). If reached, the loop is exited with exit code "stationnary" (Those points should belong to Mandelbrot set interior). A typical value is 1.e-3 calc_d2zndc2: bool If True, activates the additional computation of *d2zndz2*, needed only for the alternative normal map shading 'Milnor'. calc_orbit: bool If True, stores the value of an orbit point @ n of exit - backshift backshift: (> 0) The number of iteration backward for the stored orbit starting point Notes ===== The following complex fields will be calculated: *zn* and its derivatives (*dzndz*, *dzndc*). If calc_orbit is activated, *zn_orbit* will also be stored Exit codes are *max_iter*, *divergence*, *stationnary*. """ complex_codes = ["zn", "dzndz", "dzndc"] zn = 0 dzndz = 1 dzndc = 2 # Optionnal output fields icode = 3 d2zndc2 = -1 if calc_d2zndc2: complex_codes += ["d2zndc2"] d2zndc2 = icode icode += 1 i_znorbit = -1 if calc_orbit: complex_codes += ["zn_orbit"] i_znorbit = icode icode += 1 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.complex128 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 epscv_sq = epsilon_stationnary ** 2 deg = self.exponent iterate_once = self.from_numba_cache( "iterate_once", deg, zn, dzndz, dzndc, d2zndc2, calc_d2zndc2, max_iter, Mdiv_sq, epscv_sq ) zn_iterate = self.zn_iterate def iterate(): new_args = ( calc_orbit, i_znorbit, backshift, zn, iterate_once, zn_iterate ) args, numba_impl = self._numba_iterate_cache["calc_std_div"] if new_args == args: return numba_impl numba_impl = fs.core.numba_iterate(*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 }
[docs] @fs.utils.calc_options def newton_calc(self, *, calc_name: str = "newton_calc", subset: Optional[Fractal_array] = None, known_orders: Optional[int] = None, max_order: int = 500, max_newton: int = 20, eps_newton_cv: float = 1.e-12 ): """ Newton iterations for Mandelbrot standard set interior (power N). Parameters ========== calc_name : str The string identifier for this calculation subset : Optional `fractalshades.postproc.Fractal_array` A boolean array-like, where False no calculation is performed If `None`, all points are calculated. Defaults to `None`. known_orders : None | int list If not `None`, only the integer listed of their multiples will be candidates for the cycle order, the other rwill be disregarded. If `None` all order between 1 and `max_order` will be considered. max_order : int The maximum value tested for cycle order eps_newton_cv : float A small float to qualify the convergence of the Newton iteration Usually a fraction of a view pixel. Notes ===== .. note:: The following complex fields will be calculated: *zr* A (any) point belonging to the attracting cycle *attractivity* The cycle attractivity (for a convergent cycle it is a complex of norm < 1.) *dzrdc* Derivative of *zr* *dattrdc* Derivative of *attractivity* The following integer field will be calculated: *order* The cycle order Exit codes are 0: *max_order*, 1: *order_confirmed*. References ========== .. [1] <https://mathr.co.uk/blog/2013-04-01_interior_coordinates_in_the_mandelbrot_set.html> .. [2] <https://mathr.co.uk/blog/2014-11-02_practical_interior_distance_rendering.html> """ complex_codes = [ "zr", "attractivity", "dzrdc", "dattrdc", "_partial", "_zn" ] int_codes = ["order"] stop_codes = ["max_order", "order_confirmed"] deg = self.exponent # deg_m1 = deg - 1 # deg_m2 = deg - 2 izr = 0 idzrdz = 1 # attractivity idzrdc = 2 id2zrdzdc = 3 # dattrdc i_partial = 4 i_zn = 5 iorder = 0 reason_max_order = 1 reason_order_confirmed = 2 def set_state(): def impl(instance): instance.codes = (complex_codes, int_codes, stop_codes) instance.complex_type = np.complex128 return impl def initialize(): @numba.njit def numba_init_impl(Z, U, c): Z[4] = 1.e6 # Partial, bigger than any reasonnable np.abs(zn) return numba_init_impl zn_iterate = self.zn_iterate iterate_newton_search = self.from_numba_cache( "iterate_newton_search", deg, None, None, None, None, None, None, None, None ) def iterate(): new_args = ( known_orders, max_order, max_newton, eps_newton_cv, reason_max_order, reason_order_confirmed, izr, idzrdz, idzrdc, id2zrdzdc, i_partial, i_zn, iorder, zn_iterate, iterate_newton_search ) args, numba_impl = self._numba_iterate_cache["calc_std_div"] if new_args == args: return numba_impl numba_impl = fs.core.numba_Newton(*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, deg, zn, dzndz, dzndc, d2zndc2, calc_d2zndc2, max_iter, Mdiv_sq, epscv_sq ): """ Returns the numba implementation if exists to avoid unnecessary recompilation""" cache = self._numba_cache full_key = (key, deg, zn, dzndz, dzndc, d2zndc2, calc_d2zndc2, max_iter, Mdiv_sq, epscv_sq) try: return cache[full_key] except KeyError: deg_m1 = deg - 1 deg_m2 = deg - 2 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 _zn = Z[zn] if calc_d2zndc2: zn_m2 = _zn ** deg_m2 zn_m1 = zn_m2 * _zn zn_m = zn_m1 * _zn Z[d2zndc2] = deg * ( Z[d2zndc2] * zn_m1 + deg_m1 * Z[dzndz] * Z[dzndc] * zn_m2 ) else: zn_m1 = _zn ** deg_m1 zn_m = zn_m1 * _zn Z[dzndc] = deg * Z[dzndc] * zn_m1 + 1. Z[dzndz] = deg * Z[dzndz] * zn_m1 Z[zn] = zn_m + c if n_iter == 1: Z[dzndz] = 1. if Z[zn].real ** 2 + Z[zn].imag ** 2 > Mdiv_sq: stop_reason[0] = 1 return 1 if Z[dzndz].real ** 2 + Z[dzndz].imag ** 2 < epscv_sq: stop_reason[0] = 2 return 1 return 0 elif key == "iterate_newton_search": @numba.njit def numba_impl(d2zrdzdc, dzrdc, dzrdz, zr, c): zr_m2 = zr ** deg_m2 zr_m1 = zr_m2 * zr zr_m = zr_m1 * zr d2zrdzdc = deg * ( d2zrdzdc * zr_m1 + deg_m1 * dzrdz * dzrdc * zr_m2 ) dzrdz = deg * dzrdz * zr_m1 dzrdc = deg * dzrdc * zr_m1 + 1. zr = zr_m + c return (d2zrdzdc, dzrdc, dzrdz, zr) 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)
#============================================================================== #==============================================================================
[docs] class Perturbation_mandelbrot_N(fs.PerturbationFractal):
[docs] def __init__(self, directory: str, exponent: int): """ An arbitrary precision power-N Mandelbrot Fractal. .. math:: z_0 &= 0 \\\\ z_{n+1} &= {z_{n}}^N + c 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 """ super().__init__(directory) self.exponent = exponent # Needed for serialization # Sets default values used for postprocessing (potential) self.potential_kind = "infinity" self.potential_d = exponent 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 = "zn" # GUI 'badges' self.holomorphic = True self.implements_dzndc = "user" self.implements_fieldlines = True self.implements_newton = False self.implements_Milnor = False self.implements_interior_detection = "user" self.implements_deepzoom = True # Orbit 'on the fly' calculation self.zn_iterate = Mn_iterate(exponent) # 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 x = c0.real y = c0.imag seed_prec = mpmath.mp.prec (i, partial_dict, xr_dict ) = fsFP.perturbation_mandelbrotN_FP_loop( NP_orbit.view(dtype=np.float64), xr_detect_activated, max_orbit_iter, self.exponent, M_divergence * 2, # to be sure ref exit after close points str(x).encode('utf8'), str(y).encode('utf8'), seed_prec ) 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, epsilon_stationnary: float, BLA_eps: float=1e-6, interior_detect: bool=False, calc_dzndc: bool=True, calc_orbit: bool = False, backshift: int = 0 ): """ Perturbation iterations (arbitrary precision) for Mandelbrot standard set (power N). Parameters ========== calc_name : str The string identifier for this calculation subset: Optional `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" epsilon_stationnary: float Used only if `interior_detect` parameter is set to True A small criteria (typical range 0.01 to 0.001) used to detect earlier points belonging to a minibrot, based on dzndz1 value. If reached, the loop is exited with exit code "stationnary" BLA_eps: None | float Relative error criteria for Bilinear Approximation (default: 1.e-6) if `None` BLA is not activated. interior_detect: bool If True, activates interior points early detection. This will trigger the computation of additionnal quantities (`dzndz`), so shall be activated only when a mini fills a significative part of the view. calc_dzndc: bool If True, activates the computation of an additionnal quantities (`dzndc`), used for distance estimation plots and for shading (normal map calculation) calc_orbit: bool If True, stores the value of an orbit point @ exit - backshift backshift: int (> 0) The number of iteration backward for the stored orbit starting point """ complex_codes = ["zn"] zn = 0 code_int = 0 calc_dzndz = interior_detect if calc_dzndz: code_int += 1 complex_codes += ["dzndz"] dzndz = code_int else: dzndz = -1 if calc_dzndc: code_int += 1 complex_codes += ["dzndc"] dzndc = code_int else: dzndc = -1 if calc_orbit: code_int += 1 complex_codes += ["zn_orbit"] i_znorbit = code_int else: i_znorbit = -1 # Integer int32 fields codes "U" int_codes = ["ref_cycle_iter"] # Position in ref orbit # Stop codes stop_codes = ["max_iter", "divergence", "stationnary"] reason_max_iter = 0 reason_M_divergence = 1 reason_stationnary = 2 #---------------------------------------------------------------------- # Define the functions used for BLA approximation # BLA triggered ? BLA_activated = ( (BLA_eps is not None) and (self.dx < fs.settings.newton_zoom_level) ) nexp = self.exponent dfdz = self.from_numba_cache("dfdz", nexp, zn, dzndz, dzndc) dfdc = self.from_numba_cache("dfdc", nexp, zn, dzndz, dzndc) def set_state(): def impl(instance): instance.complex_type = np.complex128 instance.potential_M = M_divergence instance.codes = (complex_codes, int_codes, stop_codes) instance.calc_dZndz = interior_detect instance.calc_dZndc = calc_dzndc # no need for BLA_activated instance.dfdz = dfdz instance.dfdc = dfdc return impl #---------------------------------------------------------------------- # Defines initialize - jitted implementation def initialize(): new_args = (zn, dzndc, dzndz) args, numba_impl = self._numba_initialize_cache if new_args == args: return numba_impl numba_impl = fs.perturbation.numba_initialize(*new_args) self._numba_initialize_cache = (new_args, numba_impl) return numba_impl #---------------------------------------------------------------------- # Defines iterate - jitted implementation M_divergence_sq = M_divergence ** 2 epsilon_stationnary_sq = epsilon_stationnary ** 2 # Xr triggered for ultra-deep zoom xr_detect_activated = self.xr_detect_activated p_iter_zn = self.from_numba_cache( "p_iter_zn",nexp, zn, dzndz, dzndc ) p_iter_dzndz = self.from_numba_cache( "p_iter_dzndz", nexp, zn, dzndz, dzndc ) p_iter_dzndc = self.from_numba_cache( "p_iter_dzndc", nexp, zn, dzndz, dzndc ) zn_iterate = self.zn_iterate def iterate(): new_args = ( max_iter, M_divergence_sq, epsilon_stationnary_sq, reason_max_iter, reason_M_divergence, reason_stationnary, xr_detect_activated, BLA_activated, zn, dzndc, dzndz, p_iter_zn, p_iter_dzndz, p_iter_dzndc, calc_dzndc, calc_dzndz, calc_orbit, i_znorbit, backshift, zn_iterate ) args, numba_impl = self._numba_iterate_cache if new_args == args: return numba_impl numba_impl = fs.perturbation.numba_iterate(*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, nexp, zn, dzndz, dzndc): """ Returns the numba implementation if exists to avoid unnecessary recompilation""" cache = self._numba_cache full_key = (key, nexp, zn, dzndz, dzndc) try: return cache[full_key] except KeyError: # We will need also Binomial coefficients - defined as float 64 due # to numba Xrange implementation of operations C_binom = np.empty((nexp + 1), dtype=np.float64) for k in range(nexp + 1): C_binom[k] = math.comb(nexp, k) float_nexp = float(self.exponent) if key == "dfdz": @numba.njit def numba_impl(z): tmp = z # ensure same datatype as z (incl. xrange) for k in range(2, nexp): tmp = tmp * z # Has the power k return tmp * float_nexp elif key == "dfdc": @numba.njit def numba_impl(z): return 1. elif key == "p_iter_zn": @numba.njit def numba_impl(Z, ref_zn, c): # Uses a full binomial expansion - could be optimized tmp = Z[zn] * (Z[zn] + C_binom[1] * ref_zn) ref_zn_pk = ref_zn # The iterative ref_zn ** k for k in range(2, nexp): # k is the power of ref_zn # Full binomial expansion # decreasing Z[zn] power, increasing ref_zn ref_zn_pk = ref_zn_pk * ref_zn # ref_zn ** k tmp = Z[zn] * (tmp + C_binom[k] * ref_zn_pk) Z[zn] = tmp + c elif key == "p_iter_dzndz": @numba.njit def numba_impl(Z, ref_zn, ref_dzndz): # Use a full binomial expansion - could be optimized mul = (Z[zn] + C_binom[1] * ref_zn) tmp = Z[zn] * mul dtmpdz = ( Z[dzndz] * mul + Z[zn] * (Z[dzndz] + C_binom[1] * ref_dzndz) ) ref_zn_pk = ref_zn ref_dzn_pkdz = ref_dzndz for k in range(2, nexp): # Full binomial expansion ref_dzn_pkdz = ( k * ref_zn_pk * ref_dzndz ) ref_zn_pk = ref_zn_pk * ref_zn mul = (tmp + C_binom[k] * ref_zn_pk) dtmpdz = ( Z[dzndz] * mul + Z[zn] * (dtmpdz + C_binom[k] * ref_dzn_pkdz) ) tmp = Z[zn] * mul Z[dzndz] = dtmpdz elif key == "p_iter_dzndc": @numba.njit def numba_impl(Z, ref_zn, ref_dzndc): # Use a full binomial expansion - could be optimized mul = (Z[zn] + C_binom[1] * ref_zn) tmp = Z[zn] * mul dtmpdc = ( Z[dzndc] * mul + Z[zn] * (Z[dzndc] + C_binom[1] * ref_dzndc) ) ref_zn_pk = ref_zn ref_dzn_pkdc = ref_dzndc for k in range(2, nexp): # Full binomial expansion ref_dzn_pkdc = ( k * ref_zn_pk * ref_dzndc ) ref_zn_pk = ref_zn_pk * ref_zn mul = (tmp + C_binom[k] * ref_zn_pk) dtmpdc = ( Z[dzndc] * mul + Z[zn] * (dtmpdc + C_binom[k] * ref_dzn_pkdc) ) tmp = Z[zn] * mul Z[dzndc] = dtmpdc 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""" print("*** in Mandelbrot n, _ball_method", self.exponent) x = c.real y = c.imag seed_prec = mpmath.mp.prec order = fsFP.perturbation_mandelbrotN_ball_method( str(x).encode('utf8'), str(y).encode('utf8'), seed_prec, str(px).encode('utf8'), self.exponent, maxiter, M_divergence ) if order == -1: return None return order def find_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 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, " "Use 'find_any_nucleus'" ) 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)) is_ok, val = fsFP.perturbation_mandelbrotN_find_any_nucleus( str(x).encode('utf8'), str(y).encode('utf8'), seed_prec, self.exponent, order, max_newton, str(eps_cv).encode('utf8'), str(eps_pixel).encode('utf8'), ) 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 nucleus_size = fsFP.perturbation_mandelbrotN_nucleus_size_estimate( str(x).encode('utf8'), str(y).encode('utf8'), seed_prec, self.exponent, order ) nucleus_size = np.abs(nucleus_size) expo = self.exponent julia_exp = ((expo + 1.) * (expo - 1.) / expo / expo) julia_size = np.power(nucleus_size, julia_exp) return nucleus_size, julia_size #============================================================================== # 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): return super().newton_search(x, y, pix, dps, maxiter, radius_pixels)