TheAlgorithms-Python

1
"""Author Alexandre De Zotti
2

3
Draws Julia sets of quadratic polynomials and exponential maps.
4
 More specifically, this iterates the function a fixed number of times
5
 then plots whether the absolute value of the last iterate is greater than
6
 a fixed threshold (named "escape radius"). For the exponential map this is not
7
 really an escape radius but rather a convenient way to approximate the Julia
8
 set with bounded orbits.
9

10
The examples presented here are:
11
- The Cauliflower Julia set, see e.g.
12
https://en.wikipedia.org/wiki/File:Julia_z2%2B0,25.png
13
- Other examples from https://en.wikipedia.org/wiki/Julia_set
14
- An exponential map Julia set, ambiantly homeomorphic to the examples in
15
https://www.math.univ-toulouse.fr/~cheritat/GalII/galery.html
16
 and
17
https://ddd.uab.cat/pub/pubmat/02141493v43n1/02141493v43n1p27.pdf
18

19
Remark: Some overflow runtime warnings are suppressed. This is because of the
20
 way the iteration loop is implemented, using numpy's efficient computations.
21
 Overflows and infinites are replaced after each step by a large number.
22
"""
23

24
import warnings
25
from collections.abc import Callable
26
from typing import Any
27

28
import matplotlib.pyplot as plt
29
import numpy as np
30

31
c_cauliflower = 0.25 + 0.0j
32
c_polynomial_1 = -0.4 + 0.6j
33
c_polynomial_2 = -0.1 + 0.651j
34
c_exponential = -2.0
35
nb_iterations = 56
36
window_size = 2.0
37
nb_pixels = 666
38

39

40
def eval_exponential(c_parameter: complex, z_values: np.ndarray) -> np.ndarray:
41
    """
42
    Evaluate $e^z + c$.
43
    >>> eval_exponential(0, 0)
44
    1.0
45
    >>> abs(eval_exponential(1, np.pi*1.j)) < 1e-15
46
    True
47
    >>> abs(eval_exponential(1.j, 0)-1-1.j) < 1e-15
48
    True
49
    """
50
    return np.exp(z_values) + c_parameter
51

52

53
def eval_quadratic_polynomial(c_parameter: complex, z_values: np.ndarray) -> np.ndarray:
54
    """
55
    >>> eval_quadratic_polynomial(0, 2)
56
    4
57
    >>> eval_quadratic_polynomial(-1, 1)
58
    0
59
    >>> round(eval_quadratic_polynomial(1.j, 0).imag)
60
    1
61
    >>> round(eval_quadratic_polynomial(1.j, 0).real)
62
    0
63
    """
64
    return z_values * z_values + c_parameter
65

66

67
def prepare_grid(window_size: float, nb_pixels: int) -> np.ndarray:
68
    """
69
    Create a grid of complex values of size nb_pixels*nb_pixels with real and
70
     imaginary parts ranging from -window_size to window_size (inclusive).
71
    Returns a numpy array.
72

73
    >>> prepare_grid(1,3)
74
    array([[-1.-1.j, -1.+0.j, -1.+1.j],
75
           [ 0.-1.j,  0.+0.j,  0.+1.j],
76
           [ 1.-1.j,  1.+0.j,  1.+1.j]])
77
    """
78
    x = np.linspace(-window_size, window_size, nb_pixels)
79
    x = x.reshape((nb_pixels, 1))
80
    y = np.linspace(-window_size, window_size, nb_pixels)
81
    y = y.reshape((1, nb_pixels))
82
    return x + 1.0j * y
83

84

85
def iterate_function(
86
    eval_function: Callable[[Any, np.ndarray], np.ndarray],
87
    function_params: Any,
88
    nb_iterations: int,
89
    z_0: np.ndarray,
90
    infinity: float | None = None,
91
) -> np.ndarray:
92
    """
93
    Iterate the function "eval_function" exactly nb_iterations times.
94
    The first argument of the function is a parameter which is contained in
95
    function_params. The variable z_0 is an array that contains the initial
96
    values to iterate from.
97
    This function returns the final iterates.
98

99
    >>> iterate_function(eval_quadratic_polynomial, 0, 3, np.array([0,1,2])).shape
100
    (3,)
101
    >>> np.round(iterate_function(eval_quadratic_polynomial,
102
    ... 0,
103
    ... 3,
104
    ... np.array([0,1,2]))[0])
105
    0j
106
    >>> np.round(iterate_function(eval_quadratic_polynomial,
107
    ... 0,
108
    ... 3,
109
    ... np.array([0,1,2]))[1])
110
    (1+0j)
111
    >>> np.round(iterate_function(eval_quadratic_polynomial,
112
    ... 0,
113
    ... 3,
114
    ... np.array([0,1,2]))[2])
115
    (256+0j)
116
    """
117

118
    z_n = z_0.astype("complex64")
119
    for _ in range(nb_iterations):
120
        z_n = eval_function(function_params, z_n)
121
        if infinity is not None:
122
            np.nan_to_num(z_n, copy=False, nan=infinity)
123
            z_n[abs(z_n) == np.inf] = infinity
124
    return z_n
125

126

127
def show_results(
128
    function_label: str,
129
    function_params: Any,
130
    escape_radius: float,
131
    z_final: np.ndarray,
132
) -> None:
133
    """
134
    Plots of whether the absolute value of z_final is greater than
135
    the value of escape_radius. Adds the function_label and function_params to
136
    the title.
137

138
    >>> show_results('80', 0, 1, np.array([[0,1,.5],[.4,2,1.1],[.2,1,1.3]]))
139
    """
140

141
    abs_z_final = (abs(z_final)).transpose()
142
    abs_z_final[:, :] = abs_z_final[::-1, :]
143
    plt.matshow(abs_z_final < escape_radius)
144
    plt.title(f"Julia set of ${function_label}$, $c={function_params}$")
145
    plt.show()
146

147

148
def ignore_overflow_warnings() -> None:
149
    """
150
    Ignore some overflow and invalid value warnings.
151

152
    >>> ignore_overflow_warnings()
153
    """
154
    warnings.filterwarnings(
155
        "ignore", category=RuntimeWarning, message="overflow encountered in multiply"
156
    )
157
    warnings.filterwarnings(
158
        "ignore",
159
        category=RuntimeWarning,
160
        message="invalid value encountered in multiply",
161
    )
162
    warnings.filterwarnings(
163
        "ignore", category=RuntimeWarning, message="overflow encountered in absolute"
164
    )
165
    warnings.filterwarnings(
166
        "ignore", category=RuntimeWarning, message="overflow encountered in exp"
167
    )
168

169

170
if __name__ == "__main__":
171
    z_0 = prepare_grid(window_size, nb_pixels)
172

173
    ignore_overflow_warnings()  # See file header for explanations
174

175
    nb_iterations = 24
176
    escape_radius = 2 * abs(c_cauliflower) + 1
177
    z_final = iterate_function(
178
        eval_quadratic_polynomial,
179
        c_cauliflower,
180
        nb_iterations,
181
        z_0,
182
        infinity=1.1 * escape_radius,
183
    )
184
    show_results("z^2+c", c_cauliflower, escape_radius, z_final)
185

186
    nb_iterations = 64
187
    escape_radius = 2 * abs(c_polynomial_1) + 1
188
    z_final = iterate_function(
189
        eval_quadratic_polynomial,
190
        c_polynomial_1,
191
        nb_iterations,
192
        z_0,
193
        infinity=1.1 * escape_radius,
194
    )
195
    show_results("z^2+c", c_polynomial_1, escape_radius, z_final)
196

197
    nb_iterations = 161
198
    escape_radius = 2 * abs(c_polynomial_2) + 1
199
    z_final = iterate_function(
200
        eval_quadratic_polynomial,
201
        c_polynomial_2,
202
        nb_iterations,
203
        z_0,
204
        infinity=1.1 * escape_radius,
205
    )
206
    show_results("z^2+c", c_polynomial_2, escape_radius, z_final)
207

208
    nb_iterations = 12
209
    escape_radius = 10000.0
210
    z_final = iterate_function(
211
        eval_exponential,
212
        c_exponential,
213
        nb_iterations,
214
        z_0 + 2,
215
        infinity=1.0e10,
216
    )
217
    show_results("e^z+c", c_exponential, escape_radius, z_final)
218