diff --git a/book/Maths/Mandelbrot.md b/book/Maths/Mandelbrot.md deleted file mode 100644 index c2ce6ce..0000000 --- a/book/Maths/Mandelbrot.md +++ /dev/null @@ -1,17 +0,0 @@ ---- -jupytext: - formats: md:myst - text_representation: - extension: .md - format_name: myst - format_version: 0.13 - jupytext_version: 1.10.3 -kernelspec: - display_name: Python 3 (ipykernel) - language: python - name: python3 ---- - -# Mandelbrot set - -To be written. diff --git a/book/Maths/Mandelbrot/.ipynb_checkpoints/Mandelbrot-checkpoint.md b/book/Maths/Mandelbrot/.ipynb_checkpoints/Mandelbrot-checkpoint.md new file mode 100644 index 0000000..7f81e7b --- /dev/null +++ b/book/Maths/Mandelbrot/.ipynb_checkpoints/Mandelbrot-checkpoint.md @@ -0,0 +1,98 @@ +--- +jupytext: + formats: md:myst + text_representation: + extension: .md + format_name: myst + format_version: 0.13 + jupytext_version: 1.10.3 +kernelspec: + display_name: Python 3 (ipykernel) + language: python + name: python3 +--- + +# Mandelbrot set + + + +The [Mandelbrot set](https://en.wikipedia.org/wiki/Mandelbrot_set) is one of these iconic abstract object popular outside of mathematics for its aesthetic appeal. +Along with LOL cats, YouTube is filled with psychedelic videos deep diving into the incredibly complicated fractal structure of the Mandelbrot set. +Like many other fractal objects, the Mandelbrot set is often used in [Generative Arts](https://en.wikipedia.org/wiki/Generative_art), see for instance [Olbaid-ST](https://www.deviantart.com/olbaid-st/gallery/36049300/ultrafractal). + +Despite what its name suggests, the Mandelbrot set was not discovered by [Benoît Mandelbrot](https://en.wikipedia.org/wiki/Benoit_Mandelbrot), a Polish-born French-American mathematician. +It finds its origin in the field of **complex dynamics**, first studied by [Pierre Fatou](https://en.wikipedia.org/wiki/Pierre_Fatou) and [Gaston Julia](https://en.wikipedia.org/wiki/Gaston_Julia) in the early 1900s. +It was properly defined in 1978 by [Robert W. Brooks](https://en.wikipedia.org/wiki/Robert_W._Brooks) and Peter Matelski. +They were also the firsts to publish a computer-generated visualization of the Mandelbrot set, shown below. +[Adrien Douaby](https://en.wikipedia.org/wiki/Adrien_Douady) and [John H. Hubbard](https://en.wikipedia.org/wiki/John_H._Hubbard) were the first to really study the mathematical properties of the Mandelbrot set. +They are also the ones who named the set in honor of Mandelbrot for his influential work in fractal geometry. + + + +**
The first published picture of the Mandelbrot set, by R. Brooks & P. Matelski in 1978. From [Wikipedia](https://en.wikipedia.org/wiki/Mandelbrot_set).
** + +Today, even a fairly standard laptop has enough computational power to generate visualizations of the Mandelbrot set with far more details than the 1978 image. +Mathematicians and artists alike regularly come up with new ways to render this incredible mathematical object. +One such example is shown below. + +```{code-cell} ipython3 +--- +render: + image: + align: center +tags: [remove-input] +--- +import numpy as np +import numba as nb +import matplotlib +import matplotlib.pyplot as plt +from matplotlib import colors +import time + +# --> Compute the Mandelbrot set. +@nb.jit(nopython=True) +def mandelbrot_kernel(c, maxiter, horizon): + z = 0 + for i in range(maxiter): + z = z**2 + c + + if np.abs(z) > horizon: + return z, i + + return z, 0 + +@np.vectorize +def mandelbrot(c, maxiter=1000, horizon=2): + return mandelbrot_kernel(c, maxiter, horizon) + +# --> Mesh the complex plane. +cr = np.linspace(-2.25, 0.75, 3000) +ci = np.linspace(-1.25, 1.25, 2500) +c = cr[:, None] + 1j*ci[None, :] + +# --> Computation. +maxiter, horizon = 1000, 2**20 +z, n = mandelbrot(c, maxiter, horizon) + +# --> Render the Mandelbrot set. +log_horizon = np.log2(np.log(horizon)) +with np.errstate(invalid="ignore"): + M = np.nan_to_num(n + 1 - np.log2(np.log(np.abs(z))) + log_horizon) + +fig, ax = plt.subplots(1, 1, figsize=(8, 8)) + +light = colors.LightSource(azdeg=315, altdeg=10) +M = light.shade(M.T, cmap=plt.cm.gray, vert_exag=1.5, norm=colors.PowerNorm(0.3), blend_mode="hsv") + +ax.imshow(M, extent=[cr.min(), cr.max(), ci.min(), ci.max()], interpolation="bicubic") + +ax.axis(False) + +year = time.strftime("%Y") +text = ("The Mandelbrot fractal set \n" "Rendered with matploblib %s, %s - https://matplotlib.org" % (matplotlib.__version__, year)) +ax.text(cr.min() + 0.025, ci.min() + 0.025, text, color="white", fontsize=12) +``` + +```{code-cell} ipython3 + +``` diff --git a/book/Maths/Mandelbrot/.ipynb_checkpoints/binary_mandelbrot-checkpoint.md b/book/Maths/Mandelbrot/.ipynb_checkpoints/binary_mandelbrot-checkpoint.md new file mode 100644 index 0000000..5126193 --- /dev/null +++ b/book/Maths/Mandelbrot/.ipynb_checkpoints/binary_mandelbrot-checkpoint.md @@ -0,0 +1,347 @@ +--- +jupytext: + formats: md:myst + text_representation: + extension: .md + format_name: myst + format_version: 0.13 + jupytext_version: 1.10.3 +kernelspec: + display_name: Python 3 (ipykernel) + language: python + name: python3 +--- + +# Simple black and white visualization + +Before getting into advanced visualizations such as [Buddhabrot](https://en.wikipedia.org/wiki/Buddhabrot), let's start simple. +Our goal for the moment will be to produce a black and white image of the Mandelbrot set. +Points belonging to the set will be shown in black and those outside in white. +The code below shows the algorithm students are most likely to come up with. + +**
Algorithm 1: Nested `for` loops
** + +```{code-cell} ipython3 +:tags: [remove-cell] + +import numpy as np +import matplotlib.pyplot as plt +``` + +```{code-cell} ipython3 +def mandelbrot(c, maxiter=1000): + # --> Boolean mask of the Mandelbrot set. + # Assume by default that all points are in the set. + M = np.ones_like(c, dtype=bool) + + # --> Loop through the points in complex plane. + m, n = c.shape + + for i in range(m): + for j in range(n): + + # --> Initial condition. + z = 0 + + # --> Loop until | z | > 2. + for k in range(maxiter): + z = z**2 + c[i, j] + + # --> Check exit condition. + if np.abs(z) > 2: + M[i, j] = False # c[i, j] not in the set. + break + return M +``` + +The code is pretty simple. +It takes as input all the points used to discretize the complex plane (as a 2D array). +For each of them, it iterates the quadratic function for `maxiter` iterations. +If during this process $\vert z_k \vert > 2$, $c_{ij}$ does not belong to the Mandelbrot set. +If the iteration runs until `k = maxiter`, we'll assume for all intents and purposes that $c_{ij}$ is in the set (although it may not in reality but that may require a far larger number of iterations). + +Let's now benchmark this first implementation. +We'll assume $c \in \left[-2.25, 0.75\right] \times \left[ -1.25i, 1.25i\right]$ and consider a relatively coarse discretization to begin with. +The real axis is discretized with 300 points the while imaginary one is discretized with 251 points. + +```{code-cell} ipython3 +:tags: [hide-cell] + +# --> Mesh the complex plane. +cr, ci = np.linspace(-2.25, 0.75, 300), np.linspace(-1.25, 1.25, 251) +c = cr[:, None] + 1j*ci[None, :] # Broadcasting trick. +``` + +```{code-cell} ipython3 +%%timeit +mandelbrot(c, maxiter=1000) +``` + +On my laptop, this computation takes a bit more than half a minute. +A high-resolution image of the Mandelbrot set typically uses ten times as many points in both directions compared to our low-resolution discretization. +This would incur approximately a 100 fold increase in computational time (i.e. more than one hour). +While this might be acceptable for a one-shot image, deep diving into the Mandelbrot set with thousands of images to generate is out of reach. + +## Excluding points from the main cardioid + +In **Algorithm 1**, we iterate the quadratic function for all the points discretizing the complex plane. +Yet, as shown in the previous section, we know analytically that some points actually belong to the Mandelbrot set. +For any given $c$, the quadratic function admits the following fixed points + +$$ +z_* = \dfrac{1 \pm \sqrt{1-4c}}{2}. +$$ + +If $\vert z_* \vert < \dfrac{1}{2}$, these fixed points are stable and $c$ belongs to the Mandelbrot set. +This gives a simple test to check ahead of time whether the current value $c$ is or is not in the set. +We can easily modify our code to include such a test and save some precious seconds. +This is shown below. + +**
Algorithm 2: Nested `for` loops + test for the main cardiod
** + +```{code-cell} ipython3 +def mandelbrot(c, maxiter=1000): + # --> Boolean mask of the Mandelbrot set. + # Assume by default that all points are in the set. + M = np.ones_like(c, dtype=bool) + + # --> Loop through the points in complex plane. + m, n = c.shape + + for i in range(m): + for j in range(n): + + zs = 1 - np.sqrt(1 - 4*c[i, j]) + # --> Check if point is the main cardiod. + if np.abs(zs) <= 1 or np.abs(zs.conj()) <= 1: + continue + + # --> Initial condition. + z = 0 + + # --> Loop until | z | > 2. + for k in range(maxiter): + z = z**2 + c[i, j] + + # --> Check exit condition. + if np.abs(z) > 2: + M[i, j] = False # c[i, j] not in the set. + break + return M +``` + +```{code-cell} ipython3 +%%timeit +mandelbrot(c, maxiter=1000) +``` + +The same computation now takes less than 10 seconds, a 4x speed-up compared to when we didn't check if $c$ belongs to the main cardioid. +If $c$ is not in the main cardioid, similar tests could be implemented to test if it instead belongs to some **period bulbs**. +This could also save us a few additional seconds to render the image. + +## Vectorizing the code + +Although we've achieved quite a significant speed-up, one major issue with our previous code is that it'll work only if the input `c` is a `np.array`. +Hence, we cannot use it to check if a particular value of $c$ is in the Mandelbrot set or not. +In `numpy`, a function that works both for a scalar or an array as input is known as a [universal function](https://numpy.org/doc/stable/reference/ufuncs.html), or **ufunc** in short. +`np.exp`, `np.cos` and others are all ufuncs. +A universal function is a function that operates on `np.array` in an element-by-element fashion. +It also supports broadcasting, type casting and several other standard `numpy` features. +To some extent, a ufunc can be understood as a "vectorized" wrapper for a function that takes a fixed number of specific inputs and produces a fixed number of specific outputs. + +NumPy provides a simple utility, `np.vectorize`, to create user-defined ufuncs. +In what follows, we'll thus implement a function `mandelbrot` taking as input a single value of $c$ and returning whether it belongs to the Mandelbrot set or not. +We'll then make it work for array input using the decorator `@np.vectorize`. + +**
Algorithm 3: Cleaner code with `@np.vectorize`
** + +```{code-cell} ipython3 +@np.vectorize +def mandelbrot(c, maxiter=1000): + + zs = 1 - np.sqrt(1 - 4*c) + if np.abs(zs) <= 1 or np.abs(zs.conj()) <= 1: + return True + + z = 0 + for i in range(maxiter): + z = z**2 + c + if np.abs(z) > 2: + return False + return True +``` + +As you can see, this new piece of code is now much cleaner than our previous implementation. +Note that we did not really get rid of the first two nested `for` loops. +They are somehow hidden in the magic numpy is doing to transform our function. +Yet, it is now more general as it'll accept both scalar values and arrays as input. +Let's benchmark it. + +```{code-cell} ipython3 +%%timeit +mandelbrot(c, maxiter=1000) +``` + +The code now takes approximately 6 seconds to execute, another 1.5x speed-up! +Despite what the name `np.vectorize` suggests, our function is not truly vectorized. +It still operates on an element-wise fashion rather than on the whole array at once. +The small speed-up essentially comes from tricks used by the NumPy developers. +In comparison, the code below shows what a truly vectorized implementation looks like. +Its drawback though is that it can only operate on arrays and not scalar values. + +**
Algorithm 4: True vectorization in `numpy`
** + +```{code-cell} ipython3 +def mandelbrot(c, maxiter=1000): + # --> Boolean mask of the Mandelbrot set. + # Assume of all points are in the set by default. + M = np.ones_like(c, dtype=bool) + + # --> Vectorized initial condition. + z = np.zeros_like(c, dtype=complex) + + for i in range(maxiter): + # --> Iterate. + z[M] = z[M]**2 + c[M] + + # --> Check escape condition. + # Escaped points are no longer iterated. + M[np.abs(z) > 2] = False + + return M +``` + +```{code-cell} ipython3 +%%timeit +mandelbrot(c, maxiter=1000) +``` + +The true vectorized implementation is 6 times faster than the one obtained from `np.vectorize`. +Had we considered the test for the main cardioid, it might have been even faster even though we would have had to play a bit with boolean logic on arrays. + +## JIT implementation with `numba` + +Computing the Mandelbrot set essentially invovles nested `for` loops. +Yet, `for` loops are terribly inefficient in pure Python. +In these cases, Python's scientific computing ecosystem has a few extremely useful packages. +One of them is [**Numba**](https://numba.pydata.org/), an open source *just-in-time* compiler that can translate a subset of Python or NumPy code into fast machine code. +From Numba's webiste + +> Numba translates Python functions to optimized machine code at runtime using the industry-standard [LLVM](https://llvm.org/) compiler library. Numba-compiled numerical algorithms in Python can approach the speeds of C or Fortran. + +What Numba does basically is to take your naïve Python implementation and generates on-the-fly specialized machine code for different array data types and layouts to optimize performance. + +```{admonition} Be careful though... +:class: danger +At first, Numba may seem kind of magic and you may be tempted to use it on all your different pieces of code. +It may not however bring the speed-up you expected it would! +For more details about when to use or not Numba JIT capabilities, please read [Numba's documentation](https://numba.readthedocs.io/en/stable/user/5minguide.html). +``` + +Despite what we've just said, our function `mandelbrot` is luckily an excellent candidate to illustrate the massive speed-up Numba can bring. +After having imported `numba`, using its JIT capabilities is as simple as adding the decorator `@numba.jit()` to your function. +To create a **ufunc**, you can use the `@numba.vectorize` decorator. +This is illustrated below. + +**
Algorithm 5: Universal function with `numba.vectorize`
** + +```{code-cell} ipython3 +import numba + +@numba.jit(nopython=True) +def in_main_cardioid(c): + z = 1 - np.sqrt(1-4*c) + return np.abs(z) <= 1 or np.abs(np.conj(z)) <= 1 + +@numba.vectorize(nopython=True) +def mandelbrot(c, maxiter): + # --> Check if point is in main cardioid. + if in_main_cardioid(c): return True + + # --> If not, check if it is nonetheless in the + # Mandelbrot set. + z = 0 + for i in range(maxiter): + if np.abs(z) > 2: + return False + z = z**2 + c + + return True +``` + +```{code-cell} ipython3 +%%timeit +mandelbrot(c, 1000) +``` + +One limitation of ufuncs generated with `numba` is that they do not accept optional arguments (although you can hack you way into it). +You have to pass all of the arguments. +In all respect, this is a fairly some price given that our code now runs in less than 50 milliseconds, a whooping 780x speed-up compared to our original implementation! +As you increase the resolution, this speed-up will get even larger. + +**
Algorithm 5: Universal function with `numba.vectorize`
** + +```{code-cell} ipython3 +import numba + +@numba.jit(nopython=True) +def in_main_cardioid(c): + z = 1 - np.sqrt(1-4*c) + return np.abs(z) <= 1 or np.abs(np.conj(z)) <= 1 + +@numba.vectorize(nopython=True) +def mandelbrot(c, maxiter): + # --> Check if point is in main cardioid. + if in_main_cardioid(c): return True + + # --> If not, check if it is nonetheless in the + # Mandelbrot set. + x, y = c.real, c.imag + x2, y2 = x*x, y*y + for i in range(maxiter): + if (x2 + y2) > 4: + return False + + y = 2*x*y + c.imag + x = x2 - y2 + c.real + x2, y2 = x*x, y*y + + return True +``` + +```{code-cell} ipython3 +%%timeit +mandelbrot(c, 1000) +``` + +## High-resolution visualization + +We are now ready to generate our first high-resolution image of the Mandelbrot set. +We'll once again consider $c \in \left[-2.25, 0.75\right] \times \left[ -1.25i, 1.25i\right]$. +The real axis will now be discretized using 3000 points, while the imaginary one is discretized with 2501 points. +For each point (provided it is not in the main cardioid), the quadratic function will be iterated up to 1000 times to check whether the point is in the Mandelbrot set or not. +The code is shown in the cell below (you can toggle it). +The `matplotlib` code used to generate the image is fairly and won't be discussed. +Note however that we could also have used `PIL` instead, another Python package for images. + +```{code-cell} ipython3 +:tags: [hide-input] + +# --> Mesh the complex plane. +cr, ci = np.linspace(-2.25, 0.75, 3000), np.linspace(-1.25, 1.25, 2501) +c = cr[:, None] + 1j*ci[None, :] # Broadcasting trick. + +# --> Get points in/out the Mandelbrot set. +M = mandelbrot(c, 1000) + +# --> Simple binary visualization. +fig, ax = plt.subplots(1, 1, figsize=(8, 8)) + +ax.imshow( + M.T, + extent=[cr.min(), cr.max(), ci.min(), ci.max()], + cmap=plt.cm.binary) + +ax.set_aspect("equal") +ax.axis(False); +``` diff --git a/book/Maths/Mandelbrot/.ipynb_checkpoints/buddhabrot-checkpoint.md b/book/Maths/Mandelbrot/.ipynb_checkpoints/buddhabrot-checkpoint.md new file mode 100644 index 0000000..85ec803 --- /dev/null +++ b/book/Maths/Mandelbrot/.ipynb_checkpoints/buddhabrot-checkpoint.md @@ -0,0 +1,447 @@ +--- +jupytext: + formats: md:myst + text_representation: + extension: .md + format_name: myst + format_version: 0.13 + jupytext_version: 1.10.3 +kernelspec: + display_name: Python 3 (ipykernel) + language: python + name: python3 +--- + +```{code-cell} ipython3 +:tags: [remove-cell] + +import numpy as np +import numba +import matplotlib.pyplot as plt + +@numba.jit(nopython=True) +def in_main_cardioid(c): + q = (c.real - 1/4)**2 + c.imag**2 + return q*(q + (c.real - 1/4)) <= c.imag**2/4 + +@numba.jit(nopython=True) +def in_period2bulb(c): + return (c.real + 1)**2 + c.imag**2 <= 1/16 + +@numba.vectorize(nopython=True) +def mandelbrot(c, maxiter): + # --> Check if point is in main cardioid. + if in_main_cardioid(c): return True + if in_period2bulb(c): return True + + # --> If not, check if it is nonetheless in the + # Mandelbrot set. + x, y = c.real, c.imag + x2, y2 = x*x, y*y + for i in range(maxiter): + if (x2 + y2) > 4: + return False + + y = 2*x*y + c.imag + x = x2 - y2 + c.real + x2, y2 = x*x, y*y + + return True +``` + +# The Buddhabrot + +Now that we know how to make simple visualizations of the Mandelbrot set, let's turn our attention to more advanced ones. +Among all of the different ways to visualize this fractal object, the most computation-intensive one is probably the [Buddhabrot](https://en.wikipedia.org/wiki/Buddhabrot). +One of its numerous renderings is shown below. + + + +
+ +In the previous section, we made visualizations highlighting if a point $c \in \mathbb{C}$ is in the Mandelbrot (or not) or how fast the dynamical system $z_{k+1} = z^2_k + c$ escaped to infinity. +The Buddhabrot is an entirely different beast. +Rather than focusing on $c$, it depicts the probability distribution over trajectories of points escaping the Mandelbrot set, i.e. the probability distribution of the sequence $z_{k+1} = z^2_k + c$ for all $c$'s outside the Mandelbrot set. +Hence, not only do we need to determine if $c$ is the Mandelbrot set, but we also need to keep track of the whole trajectory generated by the associated quadratic function. +And we need to make sure the image looks nice! + +Rendering high resolution images of the Buddhabrot as shown above can be incredibly computation-intensive. +Millions, if not billions of unstable trajectories need to be sampled. +These computations can be even more intensive if one wants to zoom in or do an animation. +Benedikt Bitterli has a nice series of [blog posts](https://benedikt-bitterli.me/buddhabrot/) describing how he's been able to render [this video](https://www.youtube.com/watch?v=zxIcydL7wwY&ab_channel=Thunabrain). +In the rest of this section, we'll only try to produce a low-resolution image (300 x 300 pixels) of the Buddhabrot. +This is already challenging enough! +It will nonetheless show the basic techniques or ideas you can use to produce higher resolution images (caution: if you do so, the computational cost will increase drastically). + + +## Simulating a trajectory + +The Buddhabrot is the probability distribution over trajectories of points escaping the Mandelbrot. +Getting a good estimate of this probability distribution requires the computation of an extremely large number of such trajectories. +Hence, this is the portion of the code that we need to optimize as much as possible. +Fortunately, most of the work has already been done in the previous section. + +```{code-cell} ipython3 +@numba.jit(nopython=True) +def trajectory(c, maxiter): + # --> Compute trajectory. + z = np.zeros((maxiter, 2), dtype=numba.float64) + + x, y = c.real, c.imag + z[0] = x, y + x2, y2 = x*x, y*y + + for i in range(1, maxiter): + + # --> Update point. + y = 2*x*y + c.imag + x = x2 - y2 + c.real + x2, y2 = x*x, y*y + + z[i] = x, y + + # --> Check bailout condition. + if (x2 + y2) > 4: + break + + return z[:i+1] +``` + +The function `trajectory` is almost identical to the function `mandelbrot` we wrote earlier. +It takes as input the value of $c \in \mathbb{C}$ that we consider and the maximum number of iterations `maxiter` beyond which we consider $c$ to be in the Mandelbrot set. +However, it no longer returns `True` or `False`. +Instead, it returns the computed trajectory `z`. +This is pretty much the only difference with our previous piece of code. + +## Simple sampling of the trajectories + +We know how to simulate efficiently a given trajectory. +Now, how do we estimate the probability distribution? +To do so, we need to sample a large number of different $c$ values and simulate the corresponding trajectories. +If a given trajectory reaches the maximum number of iterations `maxiter`, it is rejected and a new value for $c$ is sampled. +The probability distribution of our data is then obtained using a simple 2D histogram (see `np.histogram2d` for more details). +This process is re-iterated until our estimate converges. + +The question thus boils down to **how to sample the $c$ values?** +Sampling the complex plane $\mathbb{C}$ uniformly is probably the simplest thing to do. +Yet, it would be extremely slow. +However, we do know a few things about the Mandelbrot set +- If $c$ is in the main cardioid, it is in the Mandelbrot set. +Hence, we do not even need to simulate the corresponding trajectory. +- Similarly, if $c$ is the circle centered around $x + i y$ and of radius $z$, it is in the Mandelbrot set. +Once again, we don't need to simulate the corresponding trajectory. +- If $\vert c \vert > 2$, $c$ is the not the set. +Moreover, the sequence $z_{k+1} = z_k^2 + c$ grows exponentially fast. Hence, the corresponding trajectory would not contribute much to the overall probability distribution (at least not in the region of the complex plane we are interested in). +Once again, we wouldn't need to simulate it. + +Our first sampling strategy is thus fairly simple: +1. Sample $c$ uniformly in a disk of radius 2 centered at the origin. +2. If $c$ is in the cardioid or in the period-2 bulb, reject it and go back to 1. +Otherwise, compute the corresponding sequence $\left\{ z_k \right\}$. +3. If the sequence $\left\{ z_k \right\}$ reached the maximum number of iterations allowed, reject it and go back to 1. Otherwise, add it to our dataset. + +This process is repeated until we have collected a sufficiently large number of unstable trajectories. +Once this is done, the probability distribution is estimated using a 2D histogram. +The number of bins in each direction will determine the size of the Buddhabrot image we'll render. + +The next few cells implement the different steps of our algorithm. +The code should be pretty self-explanatory. + +```{code-cell} ipython3 +def uniform_sampling(): + r = np.random.uniform(0, 2) + theta = np.random.uniform(-np.pi, np.pi) + c = r * np.exp(1j*theta) + + while in_main_cardioid(c) or in_period2bulb(c): + return uniform_sampling() + else: + return c +``` + +```{code-cell} ipython3 +def sample_trajectory(maxiter, get_new_sample): + + rejected = True + + # --> Sample new trajectory if needed. + while rejected: + # --> Sample new c value. + c = get_new_sample() + + # --> New candidate trajectory. + z = trajectory(c, maxiter) + + # --> Check whether it escaped or not. + rejected = len(z) == maxiter + else: + return z +``` + +```{code-cell} ipython3 +def generate_data(n, maxiter, get_new_sample): + + # --> Initial trajectory. + z = sample_trajectory(maxiter, get_new_sample) + + # --> Add new trajectories to the dataset. + while len(z) < n: + k = len(z) + z = np.append(z, sample_trajectory(maxiter, get_new_sample), axis=0) + else: + return z +``` + +```{code-cell} ipython3 +def compute_image( + n, + maxiter=100, + get_new_sample=uniform_sampling, + nx=301, ny=301): + + # --> Generate data. + data = generate_data(n, maxiter, get_new_sample) + + # --> Generate image. + cr, ci = np.linspace(-2, 2, nx), np.linspace(-2, 2, ny) + H, _, _ = np.histogram2d(data[:, 0], data[:, 1], bins=[cr, ci]) + + return H +``` + +Let's now generate a 300 x 300 image of the Buddhabrot. +We'll sample enough trajectories to have at least 1 000 000 points from which to estimate the probability distribution. + +```{code-cell} ipython3 +%%time +H_naive = compute_image(nsamples := 10**6, maxiter := 1000) +``` + +On my laptop, this computation takes 5 to 8 minutes depending on what I'm doing on the side. + +```{code-cell} ipython3 +--- +render: + image: + align: center +tags: [hide-input] +--- +fig, ax = plt.subplots(1, 1, figsize=(8, 8)) +ax.imshow( + H_naive, + cmap=plt.cm.inferno, + extent=[-2, 2, -2, 2], + interpolation="bicubic" +) + +ax.set_aspect("equal") +ax.set_xlim(-1.75, 1.75) +ax.set_ylim(-1.25, 1.75) +ax.axis(False); +``` + +## A smarter sampling strategy + +### Edge detection using Sobel filter + +```{code-cell} ipython3 +from scipy.ndimage import sobel + +# --> Complex plane. +nx, ny = 3001, 2501 +cr = np.linspace(-2.25, 0.75, nx) +ci = np.linspace(-1.25, 1.25, ny) +c = cr[:, None] + 1j*ci[None, :] + +# --> Compute the Mandelbrot set. +M = mandelbrot(c, maxiter) + +# --> Compute its edges. +edges = np.abs(sobel(M, axis=0)) + np.abs(sobel(M, axis=1)) +``` + +```{code-cell} ipython3 +--- +render: + image: + align: center +tags: [remove-input] +--- +# --> Plot the Mandelbrot set. +fig, axes = plt.subplots(1, 2, figsize=(16, 8)) + +axes[0].imshow( + M, + extent=[ci.min(), ci.max(), cr.min(), cr.max()], + cmap=plt.cm.binary) + +axes[0].set_title("Mandelbrot set") + +axes[1].imshow( + edges, + extent=[ci.min(), ci.max(), cr.min(), cr.max()], + cmap=plt.cm.binary) + +axes[1].set_title("Edges of the Mandelbrot set") + +for ax in axes: + ax.set_aspect("equal") + ax.axis(False); +``` + +```{code-cell} ipython3 +# --> Points from which to sample. +cs = c[edges == True].ravel() + +# --> Scale of the perturbation to be added. +Δx, Δy = cr[1] - cr[0], ci[1]-ci[0] +``` + +```{code-cell} ipython3 +def smart_sampling(cs, Δx, Δy): + + # --> Perturbation to be added to the sampled point. + Δcr = np.random.uniform(-Δx/2, Δx/2) + Δci = np.random.uniform(-Δy/2, Δy/2) + + Δc = Δcr + 1j*Δci + + # --> Value of c to be tested. + idx = np.random.randint(len(cs)) + c = cs[idx] + Δc + + while in_main_cardioid(c) or in_period2bulb(c) or np.abs(c)>2: + return smart_sampling(cs, Δx, Δy) + else: + return c +``` + +```{code-cell} ipython3 +%%time +H_edges = compute_image(nsamples, maxiter=maxiter, get_new_sample=lambda : smart_sampling(cs, Δx, Δy)) +``` + +### Plotting the Buddhabrot + +```{code-cell} ipython3 +--- +render: + image: + align: center +tags: [hide-input] +--- +fig, axes = plt.subplots(1, 2, figsize=(16, 8), sharex=True, sharey=True) + +axes[0].imshow( + H_naive, + cmap=plt.cm.inferno, + extent=[-2, 2, -2, 2], + interpolation="bicubic" +) + +axes[0].set_title("Uniform sampling") + +axes[1].imshow( + H_edges, + cmap=plt.cm.inferno, + extent=[-2, 2, -2, 2], + interpolation="bicubic" +) + +axes[1].set_title("Sampling from the edges") + +for ax in axes: + ax.set_aspect("equal") + ax.axis(False); + ax.set_xlim(-1.75, 1.75) + ax.set_ylim(-1.25, 1.75) +``` + +```{code-cell} ipython3 +--- +render: + image: + align: center +tags: [hide-input] +--- +from matplotlib import colors +p = 1/2 + +fig, axes = plt.subplots(1, 2, figsize=(16, 8), sharex=True, sharey=True) + +axes[0].imshow( + H_naive, + cmap=plt.cm.inferno, + extent=[-2, 2, -2, 2], + norm=colors.PowerNorm(p), + interpolation="bicubic" +) + +axes[0].set_title("Uniform sampling") + +axes[1].imshow( + H_edges, + cmap=plt.cm.inferno, + extent=[-2, 2, -2, 2], + norm=colors.PowerNorm(p), + interpolation="bicubic" +) + +axes[1].set_title("Sampling from the edges") + +for ax in axes: + ax.set_aspect("equal") + ax.axis(False); + ax.set_xlim(-1.75, 1.75) + ax.set_ylim(-1.25, 1.75) +``` + +## To go further + +```{code-cell} ipython3 +from joblib import Parallel, delayed, cpu_count +def compute_image( + n, + maxiter=100, + get_new_sample=uniform_sampling, + nx=301, ny=301, + n_jobs=cpu_count()): + + # --> Generate data. + if n_jobs == 1: + data = generate_data(n, maxiter, get_new_sample) + else: + data = Parallel(n_jobs=n_jobs)(delayed(generate_data)(n//n_jobs, maxiter, get_new_sample) for i in range(n_jobs)) + data = np.concatenate(data) + + # --> Generate image. + cr, ci = np.linspace(-2, 2, nx), np.linspace(-2, 2, ny) + H, _, _ = np.histogram2d(data[:, 0], data[:, 1], bins=[cr, ci]) + + return H +``` + +```{code-cell} ipython3 +%%time +H = compute_image(10**7, maxiter=maxiter, get_new_sample=lambda : smart_sampling(cs, Δx, Δy)) +``` + +```{code-cell} ipython3 +fig, ax = plt.subplots(1, 1, figsize=(8, 8)) +ax.imshow( + H, + cmap=plt.cm.inferno, + extent=[-2, 2, -2, 2], + interpolation="bicubic", + norm=colors.PowerNorm(p) +) + +ax.set_aspect("equal") +ax.set_xlim(-1.75, 1.75) +ax.set_ylim(-1.25, 1.75) +ax.axis(False); +``` + +```{code-cell} ipython3 + +``` diff --git a/book/Maths/Mandelbrot/Mandelbrot.md b/book/Maths/Mandelbrot/Mandelbrot.md new file mode 100644 index 0000000..b3ebb53 --- /dev/null +++ b/book/Maths/Mandelbrot/Mandelbrot.md @@ -0,0 +1,100 @@ +--- +jupytext: + formats: md:myst + text_representation: + extension: .md + format_name: myst + format_version: 0.13 + jupytext_version: 1.10.3 +kernelspec: + display_name: Python 3 (ipykernel) + language: python + name: python3 +--- + +# Mandelbrot set + + + +The [Mandelbrot set](https://en.wikipedia.org/wiki/Mandelbrot_set) is one of these iconic abstract object popular outside of mathematics for its aesthetic appeal. +Along with LOL cats, YouTube is filled with psychedelic videos deep diving into its incredibly complicated fractal structure. +Like many other fractal objects, the Mandelbrot set is often used in introductory classes to [Generative Art](https://en.wikipedia.org/wiki/Generative_art), see for instance [Olbaid-ST](https://www.deviantart.com/olbaid-st/gallery/36049300/ultrafractal). + +Despite what its name suggests, the Mandelbrot set was not discovered by [Benoît Mandelbrot](https://en.wikipedia.org/wiki/Benoit_Mandelbrot), a Polish-born French-American mathematician. +It finds its origin in the field of **complex dynamics**, first studied by [Pierre Fatou](https://en.wikipedia.org/wiki/Pierre_Fatou) and [Gaston Julia](https://en.wikipedia.org/wiki/Gaston_Julia) in the early 1900s. +It was properly defined in 1978 by [Robert W. Brooks](https://en.wikipedia.org/wiki/Robert_W._Brooks) and Peter Matelski. +They were also the firsts to publish a computer-generated visualization of the Mandelbrot set, shown below. +[Adrien Douaby](https://en.wikipedia.org/wiki/Adrien_Douady) and [John H. Hubbard](https://en.wikipedia.org/wiki/John_H._Hubbard) were the first to really study the mathematical properties of the Mandelbrot set. +They are also the ones who named the set in honor of Mandelbrot for his influential work in fractal geometry. + + + +**
The first published picture of the Mandelbrot set, by R. Brooks & P. Matelski in 1978. From [Wikipedia](https://en.wikipedia.org/wiki/Mandelbrot_set).
** + +This visualization was primarily limited by the computational power and printing technology available back then. +Today, even a fairly standard laptop has enough computational power to generate visualizations of the Mandelbrot set with far more details than the 1978 image. +Mathematicians and artists alike regularly come up with new ways to render this incredible mathematical object. +One such example is shown below. + +```{code-cell} ipython3 +--- +render: + image: + align: center +tags: [remove-input] +--- +import numpy as np +import numba as nb +import matplotlib +import matplotlib.pyplot as plt +from matplotlib import colors +import time + +# --> Compute the Mandelbrot set. +@nb.jit(nopython=True) +def mandelbrot_kernel(c, maxiter, horizon): + z = 0 + for i in range(maxiter): + z = z**2 + c + + if np.abs(z) > horizon: + return z, i + + return z, 0 + +@np.vectorize +def mandelbrot(c, maxiter=1000, horizon=2): + return mandelbrot_kernel(c, maxiter, horizon) + +# --> Mesh the complex plane. +cr = np.linspace(-2.25, 0.75, 3000) +ci = np.linspace(-1.25, 1.25, 2500) +c = cr[:, None] + 1j*ci[None, :] + +# --> Computation. +maxiter, horizon = 1000, 2**40 +z, n = mandelbrot(c, maxiter, horizon) + +# --> Render the Mandelbrot set. +log_horizon = np.log2(np.log(horizon)) +with np.errstate(invalid="ignore"): + M = np.nan_to_num(n + 1 - np.log2(np.log(np.abs(z))) + log_horizon) + +fig, ax = plt.subplots(1, 1, figsize=(8, 8)) + +light = colors.LightSource(azdeg=315, altdeg=10) +M = light.shade(M.T, cmap=plt.cm.gray, vert_exag=1.5, norm=colors.PowerNorm(0.3), blend_mode="hsv") + +ax.imshow(M, extent=[cr.min(), cr.max(), ci.min(), ci.max()], interpolation="bicubic") + +ax.axis(False) + +year = time.strftime("%Y") +text = ("The Mandelbrot fractal set \n" "Rendered with matploblib %s, %s." % (matplotlib.__version__, year)) +ax.text(cr.min() + 0.025, ci.min() + 0.025, text, color="white", fontsize=12); +``` + +In this chapter, we'll see how to produce different visualizations of the Mandelbrot set, from a simple black and white image all the way to the [Buddhabrot](https://en.wikipedia.org/wiki/Buddhabrot). +The algorithm itself is fairly simple and cannot really be considered part of scientific computing. +Yet, we'll see how to write clean code using `np.vectorize` and how to make it incredibly fast using `numba` just-in-time compiling capabilities. +Finally, we'll explore how to use `matplotlib` to do shaded and power normalized rendering. diff --git a/book/Maths/Mandelbrot/binary_mandelbrot.md b/book/Maths/Mandelbrot/binary_mandelbrot.md new file mode 100644 index 0000000..5330273 --- /dev/null +++ b/book/Maths/Mandelbrot/binary_mandelbrot.md @@ -0,0 +1,464 @@ +--- +jupytext: + formats: md:myst + text_representation: + extension: .md + format_name: myst + format_version: 0.13 + jupytext_version: 1.10.3 +kernelspec: + display_name: Python 3 (ipykernel) + language: python + name: python3 +--- + +# Simple visualizations + +Before getting into advanced visualizations such as [Buddhabrot](https://en.wikipedia.org/wiki/Buddhabrot), let's start simple. +Our goal for the moment will be to produce a black and white image of the Mandelbrot set. +Points belonging to the set will be shown in black and those outside in white. +The code below shows the algorithm students are most likely to come up with. + +**
Algorithm 1: Nested `for` loops
** + +```{code-cell} ipython3 +:tags: [remove-cell] + +import numpy as np +import matplotlib.pyplot as plt +``` + +```{code-cell} ipython3 +def mandelbrot(c, maxiter=1000): + # --> Boolean mask of the Mandelbrot set. + # Assume by default that all points are in the set. + M = np.ones_like(c, dtype=bool) + + # --> Loop through the points in complex plane. + m, n = c.shape + + for i in range(m): + for j in range(n): + + # --> Initial condition. + z = 0 + + # --> Loop until | z | > 2. + for k in range(maxiter): + z = z**2 + c[i, j] + + # --> Check exit condition. + if np.abs(z) > 2: + M[i, j] = False # c[i, j] not in the set. + break + return M +``` + +The code is pretty simple. +It takes as input all the points used to discretize the complex plane (as a 2D array). +For each of them, it iterates the quadratic function for `maxiter` iterations. +If during this process $\vert z_k \vert > 2$, $c_{ij}$ does not belong to the Mandelbrot set. +If the iteration runs until `k = maxiter`, we'll assume for all intents and purposes that $c_{ij}$ is in the set (although it may not in reality but that may require a far larger number of iterations). + +Let's now benchmark this first implementation. +We'll assume $c \in \left[-2.25, 0.75\right] \times \left[ -1.25i, 1.25i\right]$ and consider a relatively coarse discretization to begin with. +The real axis is discretized with 300 points the while imaginary one is discretized with 251 points. + +```{code-cell} ipython3 +:tags: [hide-cell] + +# --> Mesh the complex plane. +cr, ci = np.linspace(-2.25, 0.75, 300), np.linspace(-1.25, 1.25, 251) +c = cr[:, None] + 1j*ci[None, :] # Broadcasting trick. +``` + +```{code-cell} ipython3 +%%timeit +mandelbrot(c, maxiter=1000) +``` + +On my laptop, this computation takes a bit more than half a minute. +A high-resolution image of the Mandelbrot set typically uses ten times as many points in both directions compared to our low-resolution discretization. +This would incur approximately a 100 fold increase in computational time (i.e. more than one hour). +While this might be acceptable for a one-shot image, deep diving into the Mandelbrot set with thousands of images to generate is out of reach. + +## Excluding points from the main cardioid + +In **Algorithm 1**, we iterate the quadratic function for all the points discretizing the complex plane. +Yet, as shown in the previous section, we know analytically that some points actually belong to the Mandelbrot set. +In particular, we have analytical expressions to test whether $c$ is the main cardioid or in the period-2 bulb, both being parts of the Mandelbrot set. + +**
Algorithm 2: Nested `for` loops + test for the main cardiod
** + +```{code-cell} ipython3 +def in_main_cardioid(c): + q = (c.real - 1/4)**2 + c.imag**2 + return q*(q + (c.real -1/4)) <= c.imag**2 / 4 + +def in_period2bulb(c): + return (c.real + 1)**2 + c.imag**2 <= 1/16 + +def mandelbrot(c, maxiter=1000): + # --> Boolean mask of the Mandelbrot set. + # Assume by default that all points are in the set. + M = np.ones_like(c, dtype=bool) + + # --> Loop through the points in complex plane. + m, n = c.shape + + for i in range(m): + for j in range(n): + + # --> Check if point is the main cardioid. + if in_main_cardioid(c[i, j]): continue + if in_period2bulb(c[i, j]): continue + + # --> Initial condition. + z = 0 + + # --> Loop until | z | > 2. + for k in range(maxiter): + z = z**2 + c[i, j] + + # --> Check exit condition. + if np.abs(z) > 2: + M[i, j] = False # c[i, j] not in the set. + break + return M +``` + +```{code-cell} ipython3 +%%timeit +mandelbrot(c, maxiter=1000) +``` + +The same computation now takes less than 5 seconds, a 7x speed-up compared to when we didn't check if $c$ belongs to the main cardioid. +If $c$ is not in the main cardioid, similar tests could be implemented to test if it instead belongs to some **period bulbs**. +This could also save us a few additional seconds to render the image. + +## Vectorizing the code + +Although we've achieved quite a significant speed-up, one major issue with our previous code is that it'll work only if the input `c` is a `np.array`. +Hence, we cannot use it to check if a particular value of $c$ is in the Mandelbrot set or not. +In `numpy`, a function that works both for a scalar or an array as input is known as a [universal function](https://numpy.org/doc/stable/reference/ufuncs.html), or **ufunc** in short. +`np.exp`, `np.cos` and others are all ufuncs. +A universal function is a function that operates on `np.array` in an element-by-element fashion. +It also supports broadcasting, type casting and several other standard `numpy` features. +To some extent, a ufunc can be understood as a "vectorized" wrapper for a function that takes a fixed number of specific inputs and produces a fixed number of specific outputs. + +NumPy provides a simple utility, `np.vectorize`, to create user-defined ufuncs. +In what follows, we'll thus implement a function `mandelbrot` taking as input a single value of $c$ and returning whether it belongs to the Mandelbrot set or not. +We'll then make it work for array input using the decorator `@np.vectorize`. + +**
Algorithm 3: Cleaner code with `@np.vectorize`
** + +```{code-cell} ipython3 +@np.vectorize +def mandelbrot(c, maxiter=1000): + + if in_main_cardioid(c): return True + if in_period2bulb(c): return True + + z = 0 + for i in range(maxiter): + z = z**2 + c + if np.abs(z) > 2: + return False + return True +``` + +As you can see, this new piece of code is now much cleaner than our previous implementation. +We did not really get rid of the first two nested `for` loops though. +They are somehow hidden in the magic numpy is doing to transform our function. +Yet, it is now more general as it'll accept both scalar values and arrays as input. +Let's benchmark it. + +```{code-cell} ipython3 +%%timeit +mandelbrot(c, maxiter=1000) +``` + +The code now takes only 3 seconds to execute, another 2x speed-up! +Despite what the name `np.vectorize` suggests, our function is not truly vectorized. +It still operates on an element-wise fashion rather than on the whole array at once. +The small speed-up essentially comes from tricks used by the NumPy developers. +In comparison, the code below shows what a truly vectorized implementation looks like. +Its drawback however is that it can only operate on arrays and not scalar values. + +**
Algorithm 4: True vectorization in `numpy`
** + +```{code-cell} ipython3 +def mandelbrot(c, maxiter=1000): + # --> Boolean mask of the Mandelbrot set. + # Assume of all points are in the set by default. + M = np.ones_like(c, dtype=bool) + + # --> Vectorized initial condition. + z = np.zeros_like(c, dtype=complex) + + for i in range(maxiter): + # --> Iterate. + z[M] = z[M]**2 + c[M] + + # --> Check escape condition. + # Escaped points are no longer iterated. + M[np.abs(z) > 2] = False + + return M +``` + +```{code-cell} ipython3 +%%timeit +mandelbrot(c, maxiter=1000) +``` + +The true vectorized implementation (without the test for the cardioid and period-2 bulb) is 2 times faster than the one obtained from `np.vectorize`. +Had we considered the test for the main cardioid, it would have been even faster even though we would have had to play a bit with boolean logic on arrays. + +## JIT implementation with `numba` + +Computing the Mandelbrot set essentially involves nested `for` loops. +Yet, `for` loops are terribly inefficient in pure Python. +In these cases, Python's scientific computing ecosystem has a few extremely useful packages. +One of them is [**Numba**](https://numba.pydata.org/), an open source *just-in-time* compiler that can translate a subset of Python or NumPy code into fast machine code. +From Numba's website + +> Numba translates Python functions to optimized machine code at runtime using the industry-standard [LLVM](https://llvm.org/) compiler library. Numba-compiled numerical algorithms in Python can approach the speeds of C or Fortran. + +What Numba does basically is to take your naïve Python implementation and generates on-the-fly specialized machine code for different array data types and layouts to optimize performance. + +```{admonition} Be careful though... +:class: danger +At first, Numba may seem kind of magic and you may be tempted to use it on all your different pieces of code. +It may not however bring the speed-up you expected it would! +For more details about when to use or not Numba JIT capabilities, please read [Numba's documentation](https://numba.readthedocs.io/en/stable/user/5minguide.html). +``` + +Despite what we've just said, our function `mandelbrot` is luckily an excellent candidate to illustrate the massive speed-up Numba can bring. +After having imported `numba`, using its JIT capabilities is as simple as adding the decorator `@numba.jit()` to your function. +To create a **ufunc**, you can use the `@numba.vectorize` decorator. +This is illustrated below. + +**
Algorithm 5: Universal function with `numba.vectorize`
** + +```{code-cell} ipython3 +import numba + +@numba.jit(nopython=True) +def in_main_cardioid(c): + q = (c.real - 1/4)**2 + c.imag**2 + return q*(q + (c.real - 1/4)) <= c.imag**2/4 + +@numba.jit(nopython=True) +def in_period2bulb(c): + return (c.real + 1)**2 + c.imag**2 <= 1/16 + +@numba.vectorize(nopython=True) +def mandelbrot(c, maxiter): + # --> Check if point is in main cardioid. + if in_main_cardioid(c): return True + if in_period2bulb(c): return True + + # --> If not, check if it is nonetheless in the + # Mandelbrot set. + z = 0 + for i in range(maxiter): + if np.abs(z) > 2: + return False + z = z**2 + c + + return True +``` + +```{code-cell} ipython3 +%%timeit +mandelbrot(c, 1000) +``` + +One limitation of ufuncs generated with `numba` is that they do not accept optional arguments (although you can hack you way into it). +You have to pass all of the arguments. +In all respect, this is a fairly small price given that our code now runs in 20 milliseconds, a whooping 1700x speed-up compared to our original implementation! +As you increase the resolution, this speed-up will get even larger. + +Now that we have achieved a more than significant speed-up compared to our original implementation using fairly simple programming techniques, we can start to think about further optimizing the code. +In particular, we can get rid of complex arithmetic and transform everything into real one. +Moreover, there are some redundant multiplications we can get rid of. +The code below does just that. + +**
Algorithm 6: Further optimizations
** + +```{code-cell} ipython3 +@numba.vectorize(nopython=True) +def mandelbrot(c, maxiter): + # --> Check if point is in main cardioid. + if in_main_cardioid(c): return True + if in_period2bulb(c): return True + + # --> If not, check if it is nonetheless in the + # Mandelbrot set. + x, y = c.real, c.imag + x2, y2 = x*x, y*y + for i in range(maxiter): + if (x2 + y2) > 4: + return False + + y = 2*x*y + c.imag + x = x2 - y2 + c.real + x2, y2 = x*x, y*y + + return True +``` + +```{code-cell} ipython3 +%%timeit +mandelbrot(c, 1000) +``` + +Using this optimized code gives us an additional 2-3x speed-up, or a grand total of almost 4000x speed-up compared to the pure Python code. +As an exercise, try to figure out why writing the code this way is faster than the one provided in Algorithm 5. + +**Optimizing a code** + +When it comes to optimizing a code, particularly in terms of operations count, remember famous computer scientist [Donal Knuth](https://en.wikipedia.org/wiki/Donald_Knuth) + +> Premature optimization is the root of all evil. + +Such optimizations should really come only once you've already exhausted all the possibilities offered by simpler programming techniques. +Moreover, you should not try to optimize your whole code. +That would rapidly be a waste of time! +Instead, first come up with a fairly simple yet efficient implementation. +Then, try to identify the most important bottlenecks and optimize only these, not the rest of the code. +This will be a better use of your limited time! + +## Black and white visualization + +We are now ready to generate our first high-resolution image of the Mandelbrot set. +We'll once again consider $c \in \left[-2.25, 0.75\right] \times \left[ -1.25i, 1.25i\right]$. +The real axis will now be discretized using 3000 points, while the imaginary one is discretized with 2501 points. +For each point (provided it is not in the main cardioid), the quadratic function will be iterated up to 1000 times to check whether the point is in the Mandelbrot set or not. +The code is shown in the cell below (you can toggle it). +The `matplotlib` code used to generate the image is fairly simple and won't be discussed. +We could also have used `PIL` instead, another Python package for images. + +```{code-cell} ipython3 +--- +render: + image: + align: center +tags: [hide-input] +--- + +# --> Mesh the complex plane. +cr, ci = np.linspace(-2.25, 0.75, 3000), np.linspace(-1.25, 1.25, 2501) +c = cr[:, None] + 1j*ci[None, :] # Broadcasting trick. + +# --> Get points in/out the Mandelbrot set. +M = mandelbrot(c, 1000) + +def plot(M, extent=[-2.25, 0.75, -1.25, 1.25]): + + # --> Simple binary visualization. + fig, ax = plt.subplots(1, 1, figsize=(8, 8)) + + ax.imshow( + M.T, + extent=extent, + cmap=plt.cm.binary) + + ax.set_aspect("equal") + ax.axis(False); + +plot(M) +``` + +## Adding some colours + +I remember pretty clearly the first time I generated this black and white image of the Mandelbrot set as a student. +My code probably involved a bunch of nested `for` loops and was written in MATLAB or SciLab. +It was slow. Painfully slow. +Yet, it gave me immense joy. +YouTube was not yet a thing, and I found out only a couple of years later about all the crazy way you could use to render the Mandelbrot set. +So let's add some colours to our image! + +Until now, our piece of code only gives a yes-no answer. +Either the point is in the Mandelbrot set and the corresponding pixel is colored in black, or it is not and the pixel is colored in white. +The **escape time algorithm** is the second simplest method to render the Mandelbrot set and gives us more nuance. +Rather than simply having this yes-no answer, we'll keep track of how many iterations were needed to escape if $c$ is not part of the Mandelbrot set and color the pixel accordingly using any colormap you want. +This requires very little modification of our original code. + +**
Escape time algorithm
** + +```{code-cell} ipython3 +@numba.vectorize(nopython=True) +def mandelbrot(c, maxiter): + # --> Check if point is in main cardioid. + if in_main_cardioid(c): return maxiter + if in_period2bulb(c): return maxiter + + # --> If not, check if it is nonetheless in the + # Mandelbrot set. + x, y = c.real, c.imag + x2, y2 = x*x, y*y + for i in range(maxiter): + if (x2 + y2) > 4: + return i + + y = 2*x*y + c.imag + x = x2 - y2 + c.real + x2, y2 = x*x, y*y + + return maxiter +``` + +That's it! Literally! The escape time algorithm is as simple as replacing the `True`/`False` return with the number of iterations it took to escape. +Let's now render some Mandelbrot sets. +We'll use colormaps from the `cmasher` package. +Check out its [webpage](https://cmasher.readthedocs.io/index.html) if you want to install it. + +```{code-cell} ipython3 +--- +render: + image: + align: center +tags: [remove-input] +--- +import matplotlib.pyplot as plt +import cmasher as cmr + +# --> Mesh the complex plane. +cr, ci = np.linspace(-2.25, 0.75, 3000), np.linspace(-1.25, 1.25, 2501) +c = cr[:, None] + 1j*ci[None, :] # Broadcasting trick. + +# --> Get points in/out the Mandelbrot set. +M = mandelbrot(c, 100) + +def plot(M, extent=[-2.25, 0.75, -1.25, 1.25]): + + # --> Simple binary visualization. + fig, ax = plt.subplots(1, 1, figsize=(8, 8)) + + ax.imshow( + M.T, + extent=extent, + cmap="cmr.emergency") + + ax.set_aspect("equal") + ax.axis(False); + +plot(M) +``` + +```{code-cell} ipython3 +--- +render: + image: + align: center +tags: [remove-input] +--- +cr, ci = np.linspace(0.26-0.01, 0.26+0.21, 3000), np.linspace(0, 0.2, 3000) +c = cr[:, None] + 1j*ci[None, :] # Broadcasting trick. + +M = mandelbrot(c, 100) + +plot(M, extent=[cr.min(), cr.max(), ci.min(), ci.max()]) +``` diff --git a/book/Maths/Mandelbrot/buddhabrot.md b/book/Maths/Mandelbrot/buddhabrot.md new file mode 100644 index 0000000..63e6b93 --- /dev/null +++ b/book/Maths/Mandelbrot/buddhabrot.md @@ -0,0 +1,505 @@ +--- +jupytext: + formats: md:myst + text_representation: + extension: .md + format_name: myst + format_version: 0.13 + jupytext_version: 1.10.3 +kernelspec: + display_name: Python 3 (ipykernel) + language: python + name: python3 +--- + +```{code-cell} ipython3 +:tags: [remove-cell] + +import numpy as np +import numba +import matplotlib.pyplot as plt + +@numba.jit(nopython=True) +def in_main_cardioid(c): + q = (c.real - 1/4)**2 + c.imag**2 + return q*(q + (c.real - 1/4)) <= c.imag**2/4 + +@numba.jit(nopython=True) +def in_period2bulb(c): + return (c.real + 1)**2 + c.imag**2 <= 1/16 + +@numba.vectorize(nopython=True) +def mandelbrot(c, maxiter): + # --> Check if point is in main cardioid. + if in_main_cardioid(c): return True + if in_period2bulb(c): return True + + # --> If not, check if it is nonetheless in the + # Mandelbrot set. + x, y = c.real, c.imag + x2, y2 = x*x, y*y + for i in range(maxiter): + if (x2 + y2) > 4: + return False + + y = 2*x*y + c.imag + x = x2 - y2 + c.real + x2, y2 = x*x, y*y + + return True +``` + +# The Buddhabrot + +Now that we know how to make simple visualizations of the Mandelbrot set, let's turn our attention to more advanced ones. +Among these, the most computation-intensive one is probably the [Buddhabrot](https://en.wikipedia.org/wiki/Buddhabrot). +One of its numerous renderings is shown below. + + + +
+ +In the previous section, we made visualizations highlighting if a point $c \in \mathbb{C}$ is in the Mandelbrot (or not) or how fast the dynamical system $z_{k+1} = z^2_k + c$ escaped to infinity. +The Buddhabrot is an entirely different beast. +Rather than focusing on $c$, it depicts the probability distribution over trajectories of points escaping the Mandelbrot set, i.e. the probability distribution of the sequence $z_{k+1} = z^2_k + c$ for all $c$'s outside the Mandelbrot set. +Hence, not only do we need to determine if $c$ is the Mandelbrot set, but we also need to keep track of the whole trajectory generated by the associated quadratic function. +And we need to make sure the image looks nice as well! + +Rendering high resolution images of the Buddhabrot can be incredibly hard. +Millions, if not billions of unstable trajectories need to be sampled. +These computations are even more intensive if one wants to zoom in or create an animation. +Benedikt Bitterli has a nice series of [blog posts](https://benedikt-bitterli.me/buddhabrot/) describing how he's been able to render [this video](https://www.youtube.com/watch?v=zxIcydL7wwY&ab_channel=Thunabrain). +In the rest of this section, we'll only try to produce a low-resolution image (300 x 300 pixels) of the Buddhabrot. +This is already challenging enough! +It will nonetheless show the basic techniques or ideas you can use to produce higher resolution images (**Caution**: if you do so, you may have to let your code run for quite a while). + + +## Simulating a trajectory + +The Buddhabrot is the probability distribution over trajectories of points escaping the Mandelbrot. +Getting a good estimate of this probability distribution requires the computation of an extremely large number of such trajectories. +Hence, this is the portion of the code that we need to optimize as much as possible. +Fortunately, most of the work has already been done in the previous section. + +```{code-cell} ipython3 +@numba.jit(nopython=True) +def trajectory(c, maxiter): + # --> Compute trajectory. + z = np.zeros((maxiter, 2), dtype=numba.float64) + + x, y = c.real, c.imag + z[0] = x, y + x2, y2 = x*x, y*y + + for i in range(1, maxiter): + + # --> Update point. + y = 2*x*y + c.imag + x = x2 - y2 + c.real + x2, y2 = x*x, y*y + + z[i] = x, y + + # --> Check bailout condition. + if (x2 + y2) > 4: + break + + return z[:i+1] +``` + +The function `trajectory` is almost identical to the function `mandelbrot` we wrote earlier. +It takes as input the value of $c \in \mathbb{C}$ that we consider and the maximum number of iterations `maxiter` beyond which we consider $c$ to be in the Mandelbrot set. +However, it no longer returns `True` or `False`. +Instead, it returns the computed trajectory `z`. +This is pretty much the only difference with our previous piece of code. + +## Simple sampling of the trajectories + +We know how to simulate efficiently a given trajectory. +Now, how do we estimate the probability distribution? +To do so, we need to sample a large number of different $c$ values and simulate the corresponding trajectories. +If a given trajectory reaches the maximum number of iterations `maxiter`, it is rejected and a new value for $c$ is sampled. +The probability distribution of our data is then obtained using a simple 2D histogram (see `np.histogram2d` for more details). +This process is re-iterated until our estimate converges. + +The question thus boils down to **how to sample the $c$ values?** +Sampling the complex plane $\mathbb{C}$ uniformly is probably the simplest thing to do. +Yet, it would be extremely slow. +However, we do know a few things about the Mandelbrot set +- If $c$ is in the main cardioid, it is in the Mandelbrot set. +Hence, we do not even need to simulate the corresponding trajectory. +- Similarly, if $c$ is the period-2 bulb, it is in the Mandelbrot set. +Once again, we don't need to simulate the corresponding trajectory. +- If $\vert c \vert > 2$, $c$ is the not the set. +Moreover, the sequence $z_{k+1} = z_k^2 + c$ grows exponentially fast. The corresponding trajectory would not contribute much to the overall probability distribution. +Once again, we wouldn't need to simulate it. + +Our first sampling strategy is thus fairly simple: +1. Sample $c$ uniformly in a disk of radius 2 centered at the origin. +2. If $c$ is in the cardioid or in the period-2 bulb, reject it and go back to 1. +Otherwise, compute the corresponding sequence $\left\{ z_k \right\}$. +3. If the sequence $\left\{ z_k \right\}$ reached the maximum number of iterations allowed, reject it and go back to 1. Otherwise, add it to our dataset. + +This process is repeated until we have collected a sufficiently large number of unstable trajectories. +Once this is done, the probability distribution is estimated using a 2D histogram. +The number of bins in each direction will determine the size of the Buddhabrot image we'll render. + +The next few cells implement the different steps of our algorithm. +The code should be pretty self-explanatory. + +```{code-cell} ipython3 +def uniform_sampling(): + + # --> Uniform sampling in disk or radius 2. + r = np.random.uniform(0, 2) + theta = np.random.uniform(-np.pi, np.pi) + c = r * np.exp(1j*theta) + + # --> Accept/Reject c. + while in_main_cardioid(c) or in_period2bulb(c): + return uniform_sampling() + else: + return c +``` + +```{code-cell} ipython3 +def sample_trajectory(maxiter, get_new_sample): + + rejected = True + + # --> Sample new trajectory if needed. + while rejected: + # --> Sample new c value. + c = get_new_sample() + + # --> New candidate trajectory. + z = trajectory(c, maxiter) + + # --> Accept/Reject trajectory. + rejected = len(z) == maxiter + else: + return z +``` + +```{code-cell} ipython3 +def generate_data(n, maxiter, get_new_sample): + + # --> Initial trajectory. + z = sample_trajectory(maxiter, get_new_sample) + + # --> Add new trajectories to the dataset. + # We do so until we have at least n points + # in the dataset. + while len(z) < n: + k = len(z) + z = np.append(z, sample_trajectory(maxiter, get_new_sample), axis=0) + else: + return z +``` + +```{code-cell} ipython3 +def compute_image( + n, # Number of points to include in the dataset. + maxiter=100, # Maximum number of iterations. + get_new_sample=uniform_sampling, # How to sample c. + nx=301, ny=301): # Resolution of the final image. + + # --> Generate data. + data = generate_data(n, maxiter, get_new_sample) + + # --> Generate image. + cr, ci = np.linspace(-2, 2, nx), np.linspace(-2, 2, ny) + H, _, _ = np.histogram2d(data[:, 0], data[:, 1], bins=[cr, ci]) + + return H +``` + +Let's now generate a 300 x 300 image of the Buddhabrot. +We'll sample enough trajectories to have at least 1 000 000 points from which to estimate the probability distribution. + +```{code-cell} ipython3 +%%time +H_naive = compute_image(nsamples := 10**6, maxiter := 1000) +``` + +On my laptop, this computation takes 6 to 8 minutes depending on what I'm doing on the side. +Our estimate of the probability distribution (i.e. the "Buddhabrot") is shown below. + +```{code-cell} ipython3 +--- +render: + image: + align: center +tags: [hide-input] +--- +fig, ax = plt.subplots(1, 1, figsize=(8, 8)) +ax.imshow( + H_naive, + cmap=plt.cm.inferno, + extent=[-2, 2, -2, 2], + interpolation="bicubic" +) + +ax.set_aspect("equal") +ax.set_xlim(-1.75, 1.75) +ax.set_ylim(-1.25, 1.75) +ax.axis(False); +``` + +The Buddha starts to take shape. +Yet, the image is still pretty noisy. +They are two main reasons for that: +1. We have not sampled enough trajectories. +2. Most of the trajectories we actually sampled are relatively short ones that do not contribute much to the distribution. + +One easy fix would be to sample more trajectories, say ten millions instead of one million. +As it is, our sampling strategy is however too naïve to render the image in a reasonable amount of time. +A simple strategy would be change how we sample the $c$ values to make sure we only sample long enough trajectories. +This is what we'll explore next. + +## A smarter sampling strategy + +When redering the Mandelbrot set with the *escape time algorithm*, we observed that points close to the edges of the Mandelbrot set were the ones taking the longest time to escape. +These are points we want to sample from. +But how to do it? + +There are many ways in which we could sample these points of interest. +A simple approach would be to compute the escape time map and select only values of $c$ which take more than `k` iterations to escape. +This would be a perfectly valid strategy. +However I want to illustrate another way using an edge detection technique you can use for image processing. + +Given the black and white image of the Mandelbrot set computed in the previous section, the edges are the set of points were the surrounding pixels change from `True` to `False`. +A simple way to identify these points is to compute the gradient of the image. +To do so, one can use the [Sobel operator](https://en.wikipedia.org/wiki/Sobel_operator) already implemented as `scipy.ndimage.sobel`. +This is illustrated below. + +```{code-cell} ipython3 +from scipy.ndimage import sobel + +# --> Complex plane. +nx, ny = 3001, 2501 +cr = np.linspace(-2.25, 0.75, nx) +ci = np.linspace(-1.25, 1.25, ny) +c = cr[:, None] + 1j*ci[None, :] + +# --> Compute the Mandelbrot set. +M = mandelbrot(c, maxiter) + +# --> Compute its edges. +edges = np.abs(sobel(M, axis=0)) + np.abs(sobel(M, axis=1)) +``` + +```{code-cell} ipython3 +--- +render: + image: + align: center +tags: [remove-input] +--- +# --> Plot the Mandelbrot set. +fig, axes = plt.subplots(1, 2, figsize=(16, 8)) + +axes[0].imshow( + M, + extent=[ci.min(), ci.max(), cr.min(), cr.max()], + cmap=plt.cm.binary) + +axes[0].set_title("Mandelbrot set") + +axes[1].imshow( + edges, + extent=[ci.min(), ci.max(), cr.min(), cr.max()], + cmap=plt.cm.binary) + +axes[1].set_title("Edges of the Mandelbrot set") + +for ax in axes: + ax.set_aspect("equal") + ax.axis(False); +``` + +As you can see, convolving the black and white image of the Mandelbrot set with the Sobel filter indeed detects the edges of the set. +We can then select only these points as shown below. + +```{code-cell} ipython3 +# --> Points from which to sample. +cs = c[edges == True].ravel() + +# --> Scale of the perturbation to be added. +Δx, Δy = cr[1] - cr[0], ci[1]-ci[0] +``` + +These are the points we'll sample from. +In order to add a bit more randomness to our sampling algorithm, we'll also slightly perturb each value of $c$ we sample from this set of candidate points. +Our new sampler is implemented in the next cell. + +```{code-cell} ipython3 +def reject(c): + return in_main_cardioid(c) or in_period2bulb(c) or np.abs(c)>2 + +def smart_sampling(cs, Δx, Δy): + + # --> Perturbation to be added to the sampled point. + Δcr = np.random.uniform(-Δx/2, Δx/2) + Δci = np.random.uniform(-Δy/2, Δy/2) + + Δc = Δcr + 1j*Δci + + # --> Value of c to be tested. + idx = np.random.randint(len(cs)) + c = cs[idx] + Δc + + while reject(c): + return smart_sampling(cs, Δx, Δy) + else: + return c +``` + +Let's see how efficient our renderer becomes. +We'll once again sample trajectories until our dataset is made of one million different points. + +```{code-cell} ipython3 +%%time +H_edges = compute_image(nsamples, maxiter=maxiter, get_new_sample=lambda : smart_sampling(cs, Δx, Δy)) +``` +It now takes only a few seconds to sample the required number of points! +This is a ?x speed-up compared to our naïve sampling strategy. +We can now render the image of our estimated probability distribution. + +```{code-cell} ipython3 +--- +render: + image: + align: center +tags: [hide-input] +--- +fig, axes = plt.subplots(1, 2, figsize=(16, 8), sharex=True, sharey=True) + +axes[0].imshow( + H_naive, + cmap=plt.cm.inferno, + extent=[-2, 2, -2, 2], + interpolation="bicubic" +) + +axes[0].set_title("Uniform sampling") + +axes[1].imshow( + H_edges, + cmap=plt.cm.inferno, + extent=[-2, 2, -2, 2], + interpolation="bicubic" +) + +axes[1].set_title("Sampling from the edges") + +for ax in axes: + ax.set_aspect("equal") + ax.axis(False); + ax.set_xlim(-1.75, 1.75) + ax.set_ylim(-1.25, 1.75) +``` + +As expected, the image on the right is much cleaner than the one of the left. +It is much less noisy. +This is a direct consequence of the fact that we now sample only very good candidate trajectories which contribute significantly to the overall distribution. +This is even more visible if we apply a power-law normalization of the pixel intensities before rendering the image. +This can be done with `matplotlib.colors.PowerNorm`. + +```{code-cell} ipython3 +--- +render: + image: + align: center +tags: [hide-input] +--- +from matplotlib import colors +p = 1/2 + +fig, axes = plt.subplots(1, 2, figsize=(16, 8), sharex=True, sharey=True) + +axes[0].imshow( + H_naive, + cmap=plt.cm.inferno, + extent=[-2, 2, -2, 2], + norm=colors.PowerNorm(p), + interpolation="bicubic" +) + +axes[0].set_title("Uniform sampling") + +axes[1].imshow( + H_edges, + cmap=plt.cm.inferno, + extent=[-2, 2, -2, 2], + norm=colors.PowerNorm(p), + interpolation="bicubic" +) + +axes[1].set_title("Sampling from the edges") + +for ax in axes: + ax.set_aspect("equal") + ax.axis(False); + ax.set_xlim(-1.75, 1.75) + ax.set_ylim(-1.25, 1.75) +``` + +You now know how to render a Buddhabrot! +Even though we have massively improved the speed of our renderer by choosing a better sampling strategy, rendering a high-resolution image (e.g. 1024 x 1024 pixels) will take some time as you'll need to sample hundred of millions of trajectories before your 2D histogram converges. +Running the code over night or during a whole weekend will probably do. + +## To go further + +The code we've implemented is fairly simple and can be used as the skeleton for more advanced visualizations of the Buddhabrot. +Whether you are a student keen on taking challenges or a teacher looking for project ideas for your classes, here is a non-exhaustive list of how you can build upon this code: + +- **False colors:** The very first image of the Buddhabrot shown in this chapter is known as a [Nebulabrot](https://fr.wikipedia.org/wiki/Fichier:Nebulabrot_(5000,_500,_50).png). +It is a RGB image using techniques by NASA to render images of galaxies. +Each color channel corresponds to the probability distribution generated by trajectories with different maximum number of iterations (e.g. the red channel only considers trajectories taking less than 100 iterations to escape, the blue channel is made of trajectories taking less than 1000 iterations and the green channel by trajectories taking less than 10 000 iterations). +- **The Anti-Buddhabrot:** The Buddhabrot is the probability distribution over trajectories of points **outside** the Mandelbrot set. +Only minor modifications of the code are needed to render its dual, the Anti-Buddhabrot, i.e. the probability distribution over trajectories of points **inside** the Mandelbrot set. +An example is shown below. + + + + +
+ +- **Zooming-in:** Our sampling strategy is good enough when it comes to rendering full-scale images of the Buddhabrot. +Zooming-in is a different story. +Not all points in the vicinity of the edges of the Mandelbrot (nor in $\mathbb{C}$ for that matters) contribute equally to the probability distribution in a particular region of the complex plane. +Rendering zoomed-in images of the Buddhabrot thus necessitate a different sampling strategy. +This is nicely discussed by [Benedikt Bitterli](https://benedikt-bitterli.me/buddhabrot/). +Exploring different sampling strategies to render zoomed-in versions of the Buddhabrot would be an excellent project for anyone interested in rejection sampling, importance sampling or Markov Chain Monte-Carlo simulations! + +- **Going parallel:** Sampling different trajectories is embarrassingly parallel. +Nowadays, most computers have multiple CPUs and Python offers some utility packages such as `joblib` to exploit these. +The code below illustrate one possible way to modify the function `compute_image` to sample trajectories in parallel. + +```{code-cell} ipython3 +from joblib import Parallel, delayed, cpu_count +def compute_image( + n, + maxiter=100, + get_new_sample=uniform_sampling, + nx=301, ny=301, + n_jobs=cpu_count()): + + # --> Generate data. + if n_jobs == 1: + data = generate_data(n, maxiter, get_new_sample) + else: + data = Parallel(n_jobs=n_jobs)(delayed(generate_data)(n//n_jobs, maxiter, get_new_sample) for i in range(n_jobs)) + data = np.concatenate(data) + + # --> Generate image. + cr, ci = np.linspace(-2, 2, nx), np.linspace(-2, 2, ny) + H, _, _ = np.histogram2d(data[:, 0], data[:, 1], bins=[cr, ci]) + + return H +``` + +For more advanced project ideas, you can also think about porting some portions of the code to GPUs or combining false colors images and zooms to make a video deep-diving into the Buddhabrot! +Such a project might take some serious computation time. +Be aware that, if you want to go this way, you may need more than just your laptop... diff --git a/book/Maths/Mandelbrot/definition.md b/book/Maths/Mandelbrot/definition.md new file mode 100644 index 0000000..a260f31 --- /dev/null +++ b/book/Maths/Mandelbrot/definition.md @@ -0,0 +1,86 @@ +# Definition of the Mandelbrot set + +Let $f : \mathbb{C} \mapsto \mathbb{C}$ be the quadratic function + +$$ +f(z) = z^2 + c +$$ + +with $c \in \mathbb{C}$. +Now, consider the discrete time dynamical system generated by iterating this function with $z_0 = 0$, i.e. + +$$ +\begin{aligned} + z_0 & = 0 \\ + z_{k+1} & = z_k^2 + c \quad \forall \ k \geq 1. +\end{aligned} +$$ + +For some values of $c$, $z_k$ remains bounded as $k \to \infty$. +For others, it blows up. +This is how the Mandelbrot set is defined. +It is the set of points $c \in \mathbb{C}$ such that + +$$ +\lim_{k \to \infty} \vert z_k \vert < \infty. +$$ + +For example $c = 1$ gives rise to the sequence 0, 1, 2, 5, 25, ..., which tends to infinity. +Hence, $c = 1$ is not an element of the Mandelbrot set. +On the other hand, $c=-1$ yields 0, -1, 0, -1, 0, ..., which is bounded. +As such $c = -1$ does belong to the Mandelbrot set. + +## Some analytical results + +This definition of the Mandelbrot is not very helpful when it comes to computation as we would need to test every possible values of $c$ and iterate for long enough to see if the iteration tends to infinity or not. +Luckily, we do have a few simple analytical results which help us narrow down the portion of the complex plane we need to scan. + +First, it is easy to show by induction that if $\vert z_k \vert > 2$ for a particular $c \in \mathbb{C}$, then the iteration diverges and $c$ is not the Mandelbrot set. +This result directly implies that the Mandelbrot set is inscribed in a circle of radius $2$, i.e. any $c$ such that $\vert c \vert > 2$ is not in the set. +This already narrows down the portion of the complex plane we need to scan. + +### Main cardioid + +The most prominent piece of the Mandelbrot set is known as the *main cardioid*. +It is the region of the complex plane for which the map $f(z) = z^2 + c$ has an attracting fixed point. +It can be shown that its border is given by the set of points $c$ such that + +$$ +c = \dfrac{e^{i\theta}}{2} \left( 1 - \dfrac{e^{i\theta}}{2} \right). +$$ + +This is the equation of a cardioid. + +After some simple manipulations, we can come up with a simple test to check whether a given $c$ is in the main cardioid or not. +Let $x$ (resp. $y$) be the real (resp. imaginary) part of $c$ and $p = \sqrt{x^2 + y^2}$. +The point $c = x + i y$ is in the main cardioid (and hence in the Mandelbrot set) if + +$$ +x \leq p - 2 p^2 + \dfrac{1}{4}. +$$ + +This test can equivalently be performed without the square root. +Let $q = \left( x - \dfrac{1}{4} \right)^2 + y^2$. Then $c = x + iy$ is in the Mandelbrot set if it verifies + +$$ +q \left( q + \left(x - \dfrac{1}{4} \right) \right) \leq \dfrac{1}{4}y^2. +$$ + +This simple check will save us a tremendous amount of computation. + +### Period-2 bulb + +Another region of interest is the circle on the left of the main cardioid, attached to it as the point $c=-\dfrac{3}{4}$. +This region of the complex plane correspond to values of $c$ for which the iteration has an attracting cycle of period 2. +Using simple arguments, it can be shown that it is an actual circle of radius $\dfrac{1}{4}$ centered at $c = -1$. +As such, the following test + +$$ +\left( x + 1 \right)^2 + y^2 \leq \dfrac{1}{16} +$$ + +will let us know if $c$ is in the Mandelbrot or not. + +There are infinitely many other bulbs tangent to the main cardioid. +These are however not exact circles and we do not have such simple tests for them. +Simply testing for the main cardioid and the period-2 bulb already accounts for a significant portion of the whole Mandelbrot set and will be sufficient to massively reduce the computation time needed to render the Mandelbrot set. diff --git a/book/Physics/NBody/naive.md b/book/Physics/NBody/naive.md index 463a463..7682db4 100644 --- a/book/Physics/NBody/naive.md +++ b/book/Physics/NBody/naive.md @@ -137,7 +137,7 @@ It could also severely reduce the code readability or simply not be obvious at a In these cases, Python's scientific computing ecosystem has a few extremely useful packages. One of them is [**Numba**](https://numba.pydata.org/), an open source *just-in-time* compiler that can translate a subset of Python or NumPy code into fast machine code. -From Numba's webiste +From Numba's website > Numba translates Python functions to optimized machine code at runtime using the industry-standard [LLVM](https://llvm.org/) compiler library. Numba-compiled numerical algorithms in Python can approach the speeds of C or Fortran. diff --git a/book/Physics/NBody/simulation.md b/book/Physics/NBody/simulation.md index 1d9014e..f492818 100644 --- a/book/Physics/NBody/simulation.md +++ b/book/Physics/NBody/simulation.md @@ -179,7 +179,7 @@ HTML(anim.to_html5_video()) The code is now yours to play with! You can consider any initial arrangement of particles you'd like. -Note however that, since we implemented a **particle-particle method**, the computation rapidly gets extremely expansive as you increase the number of particles. +Note however that, since we implemented a **particle-particle method**, the computation rapidly gets extremely expensive as you increase the number of particles. You can try to see how many you can simulate in a reasonnable time on your own laptop. ## Project ideas to go further diff --git a/book/_config.yml b/book/_config.yml index ca9ce56..c839a56 100644 --- a/book/_config.yml +++ b/book/_config.yml @@ -10,7 +10,7 @@ logo: logo.png # See https://jupyterbook.org/content/execute.html execute: execute_notebooks: force - timeout: 300 + timeout: 3600 # Define the name of the latex output file for PDF builds latex: diff --git a/book/_toc.yml b/book/_toc.yml index ff78088..f4016e3 100644 --- a/book/_toc.yml +++ b/book/_toc.yml @@ -53,6 +53,11 @@ parts: chapters: - file: Maths/Overview title: Overview + - file: Maths/Mandelbrot/Mandelbrot + sections: + - file: Maths/Mandelbrot/definition + - file: Maths/Mandelbrot/binary_mandelbrot + - file: Maths/Mandelbrot/buddhabrot # - file: Maths/Mandelbrot # - file: Maths/Logistic_map # - file: Maths/Lorenz