# Mandelbrot rendering code in R

I'm trying out a new continuous histogram coloring method for Mandelbrot rendering in R, and since it requires three loops (three passes), each of them doubly nested, there should be a way to vectorize the code to speed it up.

However, I'm not yet familiar enough with vectorization in R, especially when conditionals are involved, and I would appreciate the help.

I know there's a which function which gives an array of indices of the elements which obey a condition, but I'm not sure how to make it work in this case.

Basically, I would like to make some zoom animations which is why I want the code to run faster than it does now.

If any additional improvements are possible, I would very much appreciate the advice as well.

While it's not very relevant for the question, the coloring method is as follows: instead of escape count, as usual, we are counting a kind of "cumulative convergence/divergence speed", which is defined as the sum of min(1/|z_n|^2,|z_n|^2) over the whole iterations run. It is continuous, and allows for interior coloring. The disadvantage is that we have to use floating point numbers, so very deep zooms are hard, if not impossible to render.

Here is the R code I've created, with comments:

library("jcolors")
It <- 400; #Max iterations
Es <- 10^3; #Escape radius squared (much larger than 4, intentionally)
Xm <- 1600; #x pixels
Ym <- 900; #y pixels
Mxy <- matrix(rep(1,Xm*Ym),nrow=Xm,ncol=Ym); #initiate Mandelbrot array
Hxy <- Mxy; #Histogram index array
Nb <- 40; #Number of bins
hb <- rep(0,Nb); #Histogram heights vector
Mb <- rep(0,Nb); #Bin boundaries, for later
mx <- 1:Xm;
my <- 1:Ym;
Cx <- -0.5+(-1+mx/Xm*2)*1.6; #Re(C) vector
Cy <- (-1+my/Ym*2)*0.9; #Im(C) vector
ny <- 1; #First pass loop, computing Mxy
while(ny <= Ym){nx <- 1;
while(nx <= Xm){
j <- 1;
x0 <- 0;
y0 <- 0;
s <- 1;
while(j <= It){r <- x0^2+y0^2; #|z|^2
if(r<Es){x1 <- x0^2-y0^2+Cx[nx];
y1 <- 2*x0*y0+Cy[ny]}
else{x1 <- x0; y1 <- y0; break};
x0 <- x1;
y0 <- y1;
if(r>0){s <- s+min(r,1/r)}; #this is what we measure instead of escape count
j <- j+1}
Mxy[nx,ny] <- s;
nx <- nx+1}
ny <- ny+1}
Mmax <- max(Mxy);
Mmin <- min(Mxy);
DM <- Mmax-Mmin;
dM <- DM/Nb; #Bin size
Mxy <- Mxy-Mmin; #Shifting the array values so they start with 0
ny <- 1; #Second pass loop, computing the histogram
while(ny <= Ym){nx <- 1;
while(nx <= Xm){
nb <- 1;
while(nb <= Nb){if(Mxy[nx,ny] < (nb+0.01)*dM && Mxy[nx,ny] > (nb-1.01)*dM) #Bins overlap a little to avoid missed pixels
{hb[nb] <- hb[nb]+1; #Bin height increases
Hxy[nx,ny] <- nb}; #Each pixel is assigned its bin index
nb <- nb+1}
nx <- nx+1}
ny <- ny+1}
hb <- hb/max(hb); #Normalizing the histogram
nb <- 1;  #computing the new bin boundaries
while(nb < Nb){Mb[nb+1] <- Mb[nb]+dM*hb[nb];
nb <- nb+1}
ny <- 1; #Third pass loop, scaling the array according to the histogram (multiply each interval by the bin height)
while(ny <= Ym){nx <- 1;
while(nx <= Xm){
nb <- 1;
while(nb <= Nb){
if(Hxy[nx,ny]==nb){Mxy[nx,ny] <- Mb[nb]+(Mxy[nx,ny]-(nb-1)*dM)*hb[nb]; break}
nb <- nb+1}
nx <- nx+1}
ny <- ny+1};
image(Mxy, col = jcolors_contin(palette = "rainbow", reverse=FALSE)(200), axes=FALSE, xlab="", ylab="", useRaster=TRUE)
#done


And here's the resulting image: And an example of a zoom *1.0625^(-225) at the point x=0.3602404434376143632361252444, y=-0.641313061064803174860375015. Another zoom at the same point, *1.0625^(-300) and 1700 iterations. Edit:

I accepted the answer, so here's an additional comment on the coloring method itself, for future reference.

I've read a dozen posts on Mandelbrot coloring, and usually histogram coloring is used with discrete escape count measure, which is then interpolated to form a continuous image.

Using histogram coloring with bins for continuous measure is not ideal, and leads to the loss of important details (like minibrots) and highlights less important details (the "rings" and "webs" of color you can see in the images above).

I found that the best way so far is to mix the arrays without histogram scaling (a "single bin" histogram) and with histogram scaling. Then the details look good, and the colors are rich even at very high interation counts.

Here's an example of low zoom with high iteration count It=3000 and Nb=100 bins, where I use the arithmetic mean of Mxy without histogram scaling and with it. And a 8.2E-9 zoom with the same parameters: • Are you looking for anything specific in a review or just a general review? – dfhwze Sep 4 '19 at 16:56
• @dfhwze, I would mostly like to know about vectorization, and other possible ways to increase efficiency. This is my first time on this site, but as far as I see, this kind of questions is allowed here. If I need to clarify something, I will be happy to – Yuriy S Sep 4 '19 at 16:59
• This question is fine. It's just that if you specify explicitly what you are looking for in a review, you are more likely to get that concern reviewed :) – dfhwze Sep 4 '19 at 17:01
• @dfhwze, I would like to know how to improve the speed and use R more efficiently. Because R is all about vectors and parallellization, while I'm mostly using loops here. – Yuriy S Sep 4 '19 at 17:04
• Ok I glanced over your concerns in the question. Those blockquotes are usually used to indicate the task description. – dfhwze Sep 4 '19 at 17:07

There's a couple of things you can do to improve the R-style and the efficiency of your code.

Let's first look at your while loops:

After putting your comments before the code they correspond to, removing the unnecessary semicolons, and normalising your white space, the final while loop looks like this:

# Third pass loop, scaling the array according to the histogram (multiply each interval by the bin height)
ny <- 1
while (ny <= Ym) {
nx <- 1
while (nx <= Xm) {
nb <- 1
while (nb <= Nb) {
if (Hxy[nx, ny] == nb) {
Mxy[nx, ny] <- Mb[nb] + (Mxy[nx, ny] - (nb - 1) * dM) * hb[nb]
break
}
nb <- nb + 1
}
nx <- nx + 1
}
ny <- ny + 1
}


We can rewrite this using a slightly more condensed form:

# Third pass loop, scaling the array according to the histogram (multiply each interval by the bin height)
for (ny in seq(Ym)) {
for (nx in seq(Xm)) {
for (nb in seq(Nb)) {
if (Hxy[nx, ny] == nb) {
Mxy[nx, ny] <- Mb[nb] + (Mxy[nx, ny] - (nb - 1) * dM) * hb[nb]
break
}
}
}
}


The break in there isn't necessary, it just speeds things up a bit (Hxy isn't updated during the loop, and for any given (i,j), Hxy[i,j] == nb can only be true for one value of nb). So we can remove the break. But then we can also rearrange the loop; let's put the nb loop on the outside:

# Third pass loop, scaling the array according to the histogram (multiply each
# interval by the bin height)
for (nb in seq(Nb)) {
for (ny in seq(Ym)) {
for (nx in seq(Xm)) {
if (Hxy[nx, ny] == nb) {
Mxy[nx, ny] <- Mb[nb] + (Mxy[nx, ny] - (nb - 1) * dM) * hb[nb]
}
}
}
}


There are a few ways to vectorise this. You could:

• update each row in turn

• update each column in turn

• or update the whole matrix

To update each column:

for (nb in seq(Nb)) {
for (ny in seq(Ym)) {
# find those entries that need updating
idx <- which(Hxy[, ny] == nb)
# compute the replacement value for those entries
replacement <- Mb[nb] + (Mxy[idx, ny] - (nb - 1) * dM) * hb[nb]
# update the entries
Mxy[idx, ny] <- replacement
}
}

# or alternatively:
for (nb in seq(Nb)) {
for (ny in seq(Ym)) {
Mxy[, ny] <- ifelse(
Hxy[, ny] == nb,
# value if true
Mb[nb] + (Mxy[, ny] - (nb - 1) * dM) * hb[nb],
# value if false
Mxy[, ny]
)
}
}


Importantly, a matrix is a vector. So you can actually update a whole matrix in a similar vectorised way (this uses a slightly different indexing syntax to typical matrix operations - in a 2d matrix you can index an entry by either two-coordinates (the row and column) or by one (by counting down each column in turn)):

for (nb in seq(Nb)) {
# which entries need updating:
idx <- which(Hxy == nb)
# what value should they take:
replacement <- Mb[nb] + (Mxy[idx] - (nb - 1) * dM) * hb[nb]
# update those values
Mxy[idx] <- replacement
}


[edit:] Now, we still don't have a fully vectorised version of that loop.

How could that loop over Seq(Nb) be vectorised out?

For a given value of nb we find all entries in Hxy that match nb and we refer to their position in Hxy or Mxy by idx; so whereever we use nb in the definition of replacement, we could use Hxy[idx] instead (Mb[Hxy[idx]], Hxy[idx], hb[Hxy[idx]] are all vectors, because idx is a vector); so we have:

for (nb in seq(Nb)) {
# which entries need updating:
idx <- which(Hxy == nb)
# what value should they take:
replacement <- Mb[Hxy[idx]] + (Mxy[idx] - (Hxy[idx] - 1) * dM) * hb[Hxy[idx]]
# update those values
Mxy[idx] <- replacement
}


You can now throw away those index vectors:

# begone for loop
replacement <- Mb[Hxy] + (Mxy - (Hxy - 1) * dM) * hb[Hxy]
Mxy <- replacement


You've vectorised out the whole of that final loop.

However, I suspect the main limitation on the speed of your code will be the first loop rather than the second or third; and I can't see how to readily vectorise the first loop.

You might want to split the code up into three functions, one for each of the loops, and time each of them; or use profvis to work out where your code spends most of its time (https://github.com/rstudio/profvis)

If it's of use, I modified your original code, putting the matrix computation bit into a function so I could play with it. The modified code looks like this:

mandelbrot_array <- function(
# Max iterations
n_iter,
# x pixels
n_row,
# y pixels
n_col
) {
# Escape radius squared (much larger than 4, intentionally)
Es <- 10^3
Xm <- n_row
Ym <- n_col

# initiate Mandelbrot array
Mxy <- matrix(rep(1, Xm * Ym), nrow = Xm, ncol = Ym)
# Histogram index array
Hxy <- Mxy

# Number of bins
Nb <- 40
# Histogram heights vector
hb <- rep(0, Nb)
# Bin boundaries, for later
Mb <- rep(0, Nb)
mx <- 1:Xm
my <- 1:Ym
# Re(C) vector
Cx <- -0.5 + (-1 + mx / Xm * 2) * 1.6
# Im(C) vector
Cy <- (-1 + my / Ym * 2) * 0.9

# First pass loop, computing Mxy
for (ny in seq(Ym)) {
for (nx in seq(Xm)) {
x0 <- 0
y0 <- 0
s <- 1
for (j in seq(n_iter)) {
r <- x0^2 + y0^2 #|z|^2
if (r < Es) {
x1 <- x0^2 - y0^2 + Cx[nx]
y1 <- 2 * x0 * y0 + Cy[ny]
}
else {
x1 <- x0
y1 <- y0
break
}
x0 <- x1
y0 <- y1
if (r > 0) {
s <- s + min(r, 1 / r)
} # this is what we measure instead of escape count
}
Mxy[nx, ny] <- s
}
}
Mmax <- max(Mxy)
Mmin <- min(Mxy)
DM <- Mmax - Mmin
# Bin size
dM <- DM / Nb
Mxy <- Mxy - Mmin
# Second pass loop, computing the histogram

for (nb in seq(Nb)) {
idx <- which(
Mxy < (nb + 0.01) * dM &
Mxy > (nb - 1.01) * dM
)
hb[nb] <- hb[nb] + length(idx)
Hxy[idx] <- nb
}

# Normalizing the histogram
hb <- hb / max(hb)
# computing the new bin boundaries
for (nb in seq(Nb - 1)) {
Mb[nb + 1] <- Mb[nb] + dM * hb[nb]
}
# Third pass loop, has now been vectorised completely:
Mxy <- Mb[Hxy] + (Mxy - (Hxy - 1) * dM) * hb[Hxy]

Mxy
}
$$$$

• Thank you for the improvements and the review! – Yuriy S Sep 11 '19 at 15:33
• No problem. I did look into the earlier loop in a bit more detail and may post some updated code. Although vectorised code can be used to replace some of the first loop, I couldn't get more than a three fold speed up (and had to use ~ 1000x more memory) for your actual use case. You might need to use Rcpp or something to really speed up that first loop. – Russ Hyde Sep 11 '19 at 15:46

Can we do anything about that first loop? It's currently written in a very C style, modifying single elements in the matrix Mxy. Ideally, in a vectorised implementation, you'd be able to update the whole of the matrix Mxy in one go.

Since the code is pretty complicated and slow, we might need some tests and some benchmarking to ensure that after any changes we introduce the code i) still works properly, and ii) is more efficient.

I pulled out the first loop into a function:

# Some of the variables have been renamed
compute_mxy_initial <- function(
n_iter,
n_row,
n_col,
r_escape) {
# Initial Mandelbrot array
m_xy <- matrix(1, nrow = n_row, ncol = n_col)
m_x <- seq(n_row)
m_y <- seq(n_col)

# Re(C) vector
c_x <- -0.5 + (-1 + m_x / n_row * 2) * 1.6
c_y <- (-1 + m_y / n_col * 2) * 0.9

# Compute the matrix
for (ny in seq(n_col)) {
for (nx in seq(n_row)) {
x0 <- 0
y0 <- 0
s <- 1
for (j in seq(n_iter)) {
# should it be called r2 ?
r <- x0^2 + y0^2
if (r < r_escape) {
x1 <- x0^2 - y0^2 + c_x[nx]
y1 <- 2 * x0 * y0 + c_y[ny]
} else {
# are x1 / y1 just discarded here?
x1 <- x0
y1 <- y0
break
}
x0 <- x1
y0 <- y1
if (r > 0) {
s <- s + min(r, 1 / r)
}
}
m_xy[nx, ny] <- s
}
}
m_xy
}


.. and made another function compute_mxy_test with the same body (this one will be updated).

To check that the two functions give the same results you can add a few test_that tests:

set.seed(1)

test_that(
"mxy values generated by the refactored function match the original values", {
for (i in 1:30) {
n_iter <- 1 + rpois(1, 10)
n_row <- 1 + rpois(1, 10)
n_col <- 1 + rpois(1, 10)
r_escape <- rexp(1, 1 / 1000)
expect_equal(
compute_mxy_test(n_iter, n_row, n_col, r_escape),
compute_mxy_initial(n_iter, n_row, n_col, r_escape),
tolerance = 1e-4
)
}
}
)


Anytime that you modify compute_mxy_test rerun the test. (You could alternatively, and more efficiently, precompute the expected values for the matrix and just compare the values created by compute_mxy_test to these).

You can also use the bench package to compare timings, but do this sparingly because your code is pretty slow:

timings <- bench::press(
n_iter = c(10, 100, 400),
n_row = c(10, 100, 1600),
n_col = c(10, 100, 900),
{
bench::mark(
initial = compute_mxy_initial(n_iter, n_row, n_col, r_escape = 1000),
vectorised = compute_mxy_test(n_iter, n_row, n_col, r_escape = 1000)
)
}
)


Modifications to the first loop:

What I'm trying to do here is get the for (j in seq(n_iter)) outside the positional for-loops. Unfortunately, the only way I could see to do that, was to store lots more information in RAM. You could add a data.frame that stores the coordinates, squared-radius etc for each nx,ny pair such that at each new iteration, those points with small radii get updated.

Here's the basic data structure:

compute_mxy_test <- function(
n_iter,
n_row,
n_col,
r_escape) {
# Initial Mandelbrot array
m_xy <- matrix(1, nrow = n_row, ncol = n_col)

# Re(C) vector
c_x <- -0.5 + (-1 + seq(n_row) / n_row * 2) * 1.6
c_y <- (-1 + seq(n_col) / n_col * 2) * 0.9

# construct a data-frame for use in vectorised updates
coords <- expand.grid(
nx = seq(n_row),
ny = seq(n_col)
) %>%
dplyr::mutate(
# starting values for each nx/ny pair
x0 = 0,
y0 = 0,
s = 1,
r = 0
)

# Compute the matrix

# iterate over the rows of coords and then extract nx/ny from each row
for (i in seq(nrow(coords))) {
nx <- coords[i, "nx"]
ny <- coords[i, "ny"]
x0 <- coords[i, "x0"]
y0 <- coords[i, "y0"]
s <- coords[i, "s"]
r <- coords[i, "r"]
for (j in seq(n_iter)) {
# should it be called r2 ?
if (r < r_escape) {
x1 <- x0^2 - y0^2 + c_x[nx]
y1 <- 2 * x0 * y0 + c_y[ny]
s <- ifelse(r > 0, s + min(r, 1 / r), s)
x0 <- x1
y0 <- y1
r <- x0^2 + y0^2
}
}
# update the values for x0 etc for this row
coords[i, c("x0", "y0", "s", "r")] <- c(x0, y0, s, r)
}
m_xy[cbind(coords$nx, coords$ny)] <- coords$s m_xy }  THen, still working with one row at a time, we can introduce a function that will update all the coordinate-data for that row: compute_mxy_test <- function( n_iter, n_row, n_col, r_escape) { # Initial Mandelbrot array m_xy <- matrix(1, nrow = n_row, ncol = n_col) # Re(C) vector c_x <- -0.5 + (-1 + seq(n_row) / n_row * 2) * 1.6 c_y <- (-1 + seq(n_col) / n_col * 2) * 0.9 # construct a data-frame for use in vectorised updates coords <- expand.grid( nx = seq(n_row), ny = seq(n_col) ) %>% dplyr::mutate( x0 = 0, y0 = 0, s = 1, r = 0, # note these additive offsets are used in .update_coords() cx = c_x[nx], cy = c_y[ny] ) .update_coords <- function(df) { # can assume that r >= r_escape x1 <- with(df, x0^2 - y0^2 + cx) y1 <- with(df, 2 * x0 * y0 + cy) df %>% dplyr::mutate( # pmin does pairwise min between matched elements in two equal-sized vectors s = ifelse(r > 0, s + pmin(r, 1 / r), s), x0 = x1, y0 = y1, r = x0^2 + y0^2 ) } # Compute the matrix for (i in seq(nrow(coords))) { # working with a single row, but will soon generalise my_coords <- coords[i, ] for (j in seq(n_iter)) { # only update a row if it's coords are near zero rows <- which(my_coords$r < r_escape)
my_coords[rows, ] <- .update_coords(my_coords)
}
coords[i, ] <- my_coords
}
# note the use of MATRIX[cbind(vec1, vec2)], this allows you to update
# a set of specific entries in a matrix, rather than contiguous blocks of
# entries
m_xy[cbind(coords$nx, coords$ny)] <- coords$s m_xy }  Rather than doing for(index) ... for(iteration){} and updating a single row: for a given iteration, we could determine all rows that need to be updated and then update those rows. That is we can do for (iteration){find-the-relevant-rows; update-those-rows}. The final code looks like this: compute_mxy_test <- function( n_iter, n_row, n_col, r_escape) { # Initial Mandelbrot array m_xy <- matrix(1, nrow = n_row, ncol = n_col) # Re(C) vector c_x <- -0.5 + (-1 + seq(n_row) / n_row * 2) * 1.6 c_y <- (-1 + seq(n_col) / n_col * 2) * 0.9 # construct a data-frame for use in vectorised updates coords <- expand.grid( nx = seq(n_row), ny = seq(n_col) ) %>% dplyr::mutate( x0 = 0, y0 = 0, s = 1, r = 0, cx = c_x[nx], cy = c_y[ny] ) .update_coords <- function(df) { # can assume # - nx, ny, x0, y0, s, r, cx, cy are columns in df # - that r >= r_escape x1 <- with(df, x0^2 - y0^2 + cx) y1 <- with(df, 2 * x0 * y0 + cy) df %>% dplyr::mutate( s = ifelse(r > 0, s + pmin(r, 1 / r), s), x0 = x1, y0 = y1, r = x0^2 + y0^2 ) } # Compute the matrix for (j in seq(n_iter)) { rows <- which(coords$r < r_escape)
coords[rows, ] <- .update_coords(coords[rows, ])
}

m_xy[cbind(coords$nx, coords$ny)] <- coords\$s
m_xy
}


Now, you can compare the values returned, and you can compare the speed of the original implementation, and the vectorised implementation. IN my hands, for (400, 1600, 900) the vectorised code was ~ 3x as fast as the imperative for-loop code. But it used way more memory. If you run the benchmark code above, you should be able to convince yourself of this. (But you might notice that the vectorised code is considerably slower for less demanding arguments)

Although vectorised code is (usually) faster, because it pushes more work down into a lower level, the data-structures you have to set up to make it work, and the way the code looks sometimes makes it less appealing. As I said in a comment earlier, you might be better moving this code all the way down into the C level, eg by rewriting your original for-loop implementation using Rcpp.

• This was a nice challenge, thanks. Can I reiterate that when you try to refactor code in R (or in other languages) you really need to ensure that you have tests around that code and your refactorings should be done in steps that are as small as possible. – Russ Hyde Sep 11 '19 at 18:22
• Thank you again, I'll do some tests later to see if this works for my (moderate) zooming needs. I won't pretend to understand everything you say here, because I'm not a programmer, but I understand that working with C should be better, I think most fractal renders are made in C. – Yuriy S Sep 11 '19 at 18:37
• Sorry, I didn't mean you should disregard r, you can write c code that you can call from r – Russ Hyde Sep 11 '19 at 18:40
• Oh, no problem, I got that part, R has some advantages, especially when it comes to actually rendering the pictures, and it's easy to use. Not sure if I ever manage to learn C – Yuriy S Sep 11 '19 at 18:45
• I think you'd be pleasantly surprised by how similar the code in compute_mxy_initial` is to an equivalent C++ function. Try working through a few of the Rcpp exercises in the book 'Advanced R': adv-r.had.co.nz/Rcpp.html ; perhaps someone else could do an R -> C++ refactoring for just that section. (The figures look great btw) – Russ Hyde Sep 12 '19 at 9:18