Source code for fractalshades.models.collatz

# -*- coding: utf-8 -*-
import numpy as np
import numba

import fractalshades as fs
import fractalshades.utils as fsutils


[docs] class Collatz(fs.Fractal):
[docs] def __init__(self, directory): """ The Collatz fracal: .. math:: z_0 &= c \\\\ z_{n+1} &= 0.25 (2 + 7 z_n - (2+5z_n) \cos(\pi z_n)) This fractal is linked to the Syracuse conjecture. Parameters ========== directory : str Path for the working base directory """ super().__init__(directory) self.potential_kind = "transcendent"
[docs] @fsutils.calc_options def base_calc(self, *, calc_name: str, subset, max_iter: int, M_divergence: float, epsilon_stationnary: float, ): """ Basic iterations for Mandelbrot standard set (power 2). 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" 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 Notes ===== The following complex fields will be calculated: *zn* and its derivatives (*dzndz*, *dzndc*, *d2zndc2*). Exit codes are *max_iter*, *divergence*, *stationnary*. """ complex_codes = ["zn", "dzndz"] 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 return impl def initialize(): @numba.njit def numba_init_impl(Z, U, c): Z[1] = 1. Z[0] = c return numba_init_impl def iterate(): Mdiv_sq = self.M_divergence ** 2 epscv_sq = self.epsilon_stationnary ** 2 max_iter = self.max_iter @numba.njit def numba_impl(c, Z, U, stop_reason): n_iter = 0 while True: n_iter += 1 if n_iter >= max_iter: stop_reason[0] = 0 break cos0 = np.cos(np.pi * Z[0]) sin0 = np.sin(np.pi * Z[0]) Z[1] = 0.25 * Z[1] * ( 7. - 5. * cos0 + (2. + 5. * Z[0]) * np.pi * sin0 ) Z[0] = 0.25 * (2. + 7. * Z[0] - (2. + 5. * Z[0]) * cos0) if n_iter == 1: Z[1] = 1. if Z[0].real**2 + Z[0].imag ** 2 > Mdiv_sq: stop_reason[0] = 1 break if Z[1].real**2 + Z[1].imag ** 2 < epscv_sq: stop_reason[0] = 2 break # End of while loop return n_iter return numba_impl return { "set_state": set_state, "initialize": initialize, "iterate": iterate }