# Playing With Julia (Set)

Viento, me pongo en movimiento y hago crecer las olas del mar que tienes dentro (Tercer Movimiento: Lo de Dentro, Extremoduro)

I really enjoy drawing complex numbers: it is a huge source of entertainment for me. In this experiment I play with the Julia Set, another beautiful fractal like this one. This is what I have done:

• Choosing the function `f(z)=exp(z3)-0.621`
• Generating a grid of complex numbers with both real and imaginary parts in [-2, 2]
• Iterating `f(z)` over the grid a number of times so `zn+1 = f(zn)`
• Drawing the resulting grid as I did here
• Gathering all plots into a GIF with ImageMagick as I did in my previous post: each frame corresponds to a different number of iterations

This is the result:

I love how easy is doing difficult things in R. You can play with the code changing `f(z)` as well as color palettes. Be ready to get surprised:

```library(ggplot2)
library(dplyr)
library(RColorBrewer)
dir.create("output")
setwd("output")
f = function(z,c) exp(z^3)+c
# Grid of complex
z0 <- outer(seq(-2, 2, length.out = 1200),1i*seq(-2, 2, length.out = 1200),'+') %>% c()
opt <-  theme(legend.position="none",
panel.background = element_rect(fill="white"),
plot.margin=grid::unit(c(1,1,0,0), "mm"),
panel.grid=element_blank(),
axis.ticks=element_blank(),
axis.title=element_blank(),
axis.text=element_blank())
for (i in 1:35)
{
z=z0
# i iterations of f(z)
for (k in 1:i) z <- f(z, c=-0.621)
df=data.frame(x=Re(z0),
y=Im(z0),
z=as.vector(exp(-Mod(z)))) %>% na.omit()
p=ggplot(df, aes(x=x, y=y, color=z)) +
geom_tile() +
scale_x_continuous(expand=c(0,0))+
scale_y_continuous(expand=c(0,0))+
ggsave(plot=p, file=paste0("plot", stringr::str_pad(i, 4, pad = "0"),".png"), width = 1.2, height = 1.2)
}
# Place the exact path where ImageMagick is installed
system('"C:\\Program Files\\ImageMagick-6.9.3-Q16\\convert.exe" -delay 20 -loop 0 *.png julia.gif')
# cleaning up
file.remove(list.files(pattern=".png"))
```

# Zooming

You don’t have to be beautiful to turn me on (Kiss, Prince)

I discovered recently how easy is to create GIFs with R using ImageMagick and I feel like a kid with a new toy. To begin this new era of my life as R programmer I have done this:

First of all, read this article: it explains very well how to start doing GIFs from scratch. The one I have done is inspired in this previous post where I take a set of complex numbers to transform and color it using HSV technique. In this case I use this next transformation: `f(z)= -Im(z)+(Re(z)+0.5*Im(z))*1i`

Modifying the range of `Real` and `Imaginary` parts of complex numbers I obtain the zooming  effect. The code is very simple. Play with it changing the transformation or the animation options. Send me your creations, I would love to see them:

```library(dplyr)
library(ggplot2)
dir.create("output")
setwd("output")
id=1 # label tO name plots
for (i in seq(from=320, to=20, length.out = 38)){
z=outer(seq(from = -i, to = i, length.out = 300),1i*seq(from = -i, to = i, length.out = 500),'+') %>% c()
z0=z
for (k in 1:100) z <- -Im(z)+(Re(z)+0.5*Im(z))*1i
df=data.frame(x=Re(z0),
y=Im(z0),
h=(Arg(z)<0)*1+Arg(z)/(2*pi),
s=(1+sin(2*pi*log(1+Mod(z))))/2,
v=(1+cos(2*pi*log(1+Mod(z))))/2) %>% mutate(col=hsv(h,s,v))
ggplot(df, aes(x, y)) +
geom_tile(fill=df\$col)+
scale_x_continuous(expand=c(0,0))+
scale_y_continuous(expand=c(0,0))+
labs(x=NULL, y=NULL)+
theme(legend.position="none",
panel.background = element_rect(fill="white"),
plot.margin=grid::unit(c(1,1,0,0), "mm"),
panel.grid=element_blank(),
axis.ticks=element_blank(),
axis.title=element_blank(),
axis.text=element_blank())
id=id+1
}
system('"C:\\Program Files\\ImageMagick-6.9.3-Q16\\convert.exe" -delay 10 -loop 0 -duplicate 1,-2-1 *.png zooming.gif')
# cleaning up
file.remove(list.files(pattern=".png"))
```

# butteRfly

Float like a butterfly, sting like a bee (Muhammad Ali)

The Butterfly Curve was discovered by Temple H. Fay when he was in Southern University, Mississippi, and rapidly gained the attention of students and mathematicians because of its beautiful simmetry. Small dots of this plot are generated according to parametric equations of the Butterfly Curve. Big dots are randomdly distributed over the canvas:

This is the code to create butterflies:

```library(ggplot2)
npoints=500
npointsb=1200
opt=theme(legend.position="none",
panel.background = element_blank(),
panel.grid = element_blank(),
axis.ticks=element_blank(),
axis.title=element_blank(),
axis.text =element_blank())
t=seq(0,10*pi,length=npointsb)
butterfly=data.frame(x=sin(t)*(exp(1)^cos(t)-2*cos(4*t)-(sin(t/12))^5), y=cos(t)*(exp(1)^cos(t)-2*cos(4*t)-(sin(t/12))^5), s=runif(npointsb, min=.1, max=10), f=factor(sample(1:10,npointsb,TRUE)), a=runif(npointsb,min=.1, max=.4))
points=data.frame(x=runif(npoints,-4,4), y=runif(npoints,-3,5), s=runif(npoints,min=30, max=50), f=factor(sample(1:10,npoints,TRUE)), a=runif(npoints,min=.05, max=.15))
data=rbind(butterfly, points)
ggplot(data, aes(x, y, colour=f))+geom_point(alpha=data\$a,size=data\$s)+opt
```

# The Zebra Of Riemann

Mathematics is the art of giving the same name to different things (Henri Poincare)

Many surveys among experts point that demonstration of the Riemann Hypothesis is the most important pending mathematical issue in this world. This hypothesis is related to Riemann zeta function, which is supossed to be zero only for those complex whose real part is equal to 1/2 (this is the conjecture of Riemann itself). Confirming the conjecture would imply deep consequences in prime numbers teory and also in our knowledge of their properties. Next plot represents the argument of zeta function over complex with real part between -75 and 5 (x axis) and imaginary part between -40 and 40 (y axis). An explanation of this kind of graph can be found here. Does not it remind you of something?

Solving the hypothesis of Riemann is one of the seven Clay Mathematics Institute Millennium Prize Problems. Maybe zebras keep the secret to solve it in their skins.

This is the code to plot the graph:

```require(pracma)
z=outer(seq(-75, 5, by =.1),1i*seq(-40, 40, by =.1),'+')
z=apply(z, c(1,2), function(x) Arg(zeta(x)))
opt=theme(legend.position="none",
panel.background = element_blank(),
panel.grid = element_blank(),
axis.ticks=element_blank(),
axis.title=element_blank(),
axis.text =element_blank())
z=data.frame(expand.grid(x=seq(ncol(z)), y=seq(nrow(z))), z=as.vector(z))
require(ggplot2)
ggplot(z, aes(x=x, y=y, color=z)) + geom_tile() + scale_colour_gradientn(colours=c("black","white", "gray80", "gray40")) + opt
```

# floweR

It is the time you have wasted for your rose that makes your rose so important (Antoine de Saint-Exupéry, The Little Prince)

Yesterday I found a package called `circular` and I could not suppress to do this:

```require(circular)
for (i in 1:8)
{
rose.diag(circular(runif(5000, 0, 2*pi)),
bins = (25-round(abs(jitter(0, amount=2*i)))),
axes=FALSE,
border=rgb(255,50+204*((i-1)/8),50+204*((i-1)/8), alpha=255, max=255),
ticks = FALSE,
col=rgb(250,0+204*((i-1)/8),0+204*((i-1)/8), alpha=255, max=255),
control.circle=circle.control(lty=0),
shrink=0.22+(i-1)*(2-0.22)/8)
par(new=TRUE)
}
for (i in 1:150)
{
q = runif(1)*pi*2
r = sqrt(runif(1))
x = (0.18*r)* cos(q)
y = (0.18*r)* sin(q)
points(x, y,
col = rgb(255 ,sample(80:120, 1), 0, alpha= 95, max=255),
bg = rgb(255 ,sample(100:200, 1), 0, alpha= 95, max=255), pch=21, cex=2)
}
```

# Blurry Fractals

Beauty is the first test; there is no permanent place in the world for ugly mathematics (G. H. Hardy)

Newton basin fractals are the result of iterating Newton’s method to find roots of a polynomial over the complex plane. It maybe sound a bit complicated but is actually quite simple to understand. Those who would like to read some more about Newton basin fractals can visit this page.

This fractals are very easy to generate in R and produce very nice images. Making a small number of iterations, resulting images seems to be blurred when are represented with tile geometry in ggplot. Combined with palettes provided by RColorBrewer give rise to very interesting images. Here you have some examples:

Result for `f(z)=z3-1` and palette equal to `Set3`:Result for `f(z)=z4+z-1` and palette equal to `Paired`:Result for `f(z)=z5+z3+z-1` and palette equal to `Dark2`:Here you have the code. If you generate nice pictures I will be very grateful if you send them to me:

```library(ggplot2)
library(numDeriv)
library(RColorBrewer)
library(gridExtra)
## Polynom: choose only one or try yourself
f  <- function (z) {z^3-1}        #Blurry 1
#f  <- function (z) {z^4+z-1}     #Blurry 2
#f  <- function (z) {z^5+z^3+z-1} #Blurry 3
z <- outer(seq(-2, 2, by = 0.01),1i*seq(-2, 2, by = 0.01),'+')
for (k in 1:5) z <- z-f(z)/matrix(grad(f, z), nrow=nrow(z))
## Supressing texts, titles, ticks, background and legend.
opt <- theme(legend.position="none",
panel.background = element_blank(),
axis.ticks=element_blank(),
axis.title=element_blank(),
axis.text =element_blank())
z <- data.frame(expand.grid(x=seq(ncol(z)), y=seq(nrow(z))), z=as.vector(exp(-Mod(f(z)))))
# Create plots. Choose a palette with display.brewer.all()
p1 <- ggplot(z, aes(x=x, y=y, color=z)) + geom_tile() + scale_colour_gradientn(colours=brewer.pal(8, "Paired")) + opt
p2 <- ggplot(z, aes(x=x, y=y, color=z)) + geom_tile() + scale_colour_gradientn(colours=brewer.pal(7, "Paired")) + opt
p3 <- ggplot(z, aes(x=x, y=y, color=z)) + geom_tile() + scale_colour_gradientn(colours=brewer.pal(6, "Paired")) + opt
p4 <- ggplot(z, aes(x=x, y=y, color=z)) + geom_tile() + scale_colour_gradientn(colours=brewer.pal(5, "Paired")) + opt
# Arrange four plots in a 2x2 grid
grid.arrange(p1, p2, p3, p4, ncol=2)
```

# Warholing Grace With Clara

Do not believe anything: what artists really do is to hang around all day (Paco de Lucia)

Andy Warhol was mathematician. At least, he knew how clustering algorithms work. I am pretty sure of this after doing this experiment.  First of all, let me introduce you to the breathtaking Grace Kelly:

In my previous post I worked also with images showing how simple is to operate with them since they are represented by matrices. This is another example of this. Third dimension of an image matrix is an 3D array representing color of pixels in (r, g, b) format. Applying a cluster algorithm over this information generates groups of pixels with similar color. I used `cluster` package and because of the high size of picture I decided to use clara algorithm which is extremely fast. Apart of its high speed, another advantage of clara is that clusters are represented by real elements of the population, called medoids, instead of being by average individuals as k-means do. It fits very well with my purposes because once clusters are calculated I only have to change each pixel by its medoid and plot it. Setting clara to divide pixels into 2 groups, generates a 2 colored image. Setting it to 3 groups, generates a 3 colored one and so on. Following, you can find results from 2 to 7 groups:

Working with samples can be a handicap, maybe less important than the speed it produces. Sometimes images generated by n groups seems to be worse fitted than the one generated by n-1 groups. You can see it in this video, where results from 1 to 60 groups are presented sequentially. It only takes 42 seconds.

Here you have the code. Feel free to warholing:

```library("biOps")
library("abind")
library("reshape")
library("reshape2")
library("cluster")
library("sp")
#######################################################################################################
#Initialization
#######################################################################################################
plot(x)
#######################################################################################################
#Data
#######################################################################################################
data <- merge(merge(melt(x[,,1]), melt(x[,,2]), by=c("X1", "X2")), melt(x[,,3]), by=c("X1", "X2"))
colnames(data) <- c("X1", "X2", "r", "g", "b")
#######################################################################################################
#Clustering
#######################################################################################################
colors <- 5
clarax <- clara(data[,3:5], colors)
datacl   <- data.frame(data, clarax\$cluster)
colnames(claradf2) <- c("r", "g", "b", "clarax.cluster")
colnames(claradf) <-c("clarax.cluster", "X1", "X2", "r.x", "g.x", "b.x", "r.y", "g.y", "b.y")
x1<-acast(datac[,c(2,3,7)], X1~X2, value.var="r.y")
x2<-acast(datac[,c(2,3,8)], X1~X2, value.var="g.y")
x3<-acast(datac[,c(2,3,9)], X1~X2, value.var="b.y")
warhol <- do.call(abind, c(list(x1,x2,x3), along = 3))
plot(imagedata(warhol))
writeJpeg(paste("Warholing", as.character(colors), ".jpg", sep=""), imagedata(warhol))
```

# Face To Face With Marilyn Monroe

Symmetry is what we see at a glance (Blaise Pascal)

Ladies and gentlement, the beautiful Marilyn Monroe:

There are several image processing packages in R. In this experiment I used `biOps`, which turns images into 3D matrices. The third dimension is a 3-array corresponding to (r, g, b) color of pixel defined by two other dimensions. Since images are defined by matrices, you can do simple operations to produce interesting images from the original one. For example, this is what happens swaping last half of columns by first one and preserving its order:

It is also very simple to generate two artificial symmetrical faces matching each half with itself as a mirror. By the way, I prefer the first one: maybe two sexy moles are better than just one.

Let’s introduce a bit of randomness. This is what happens when you take an uniform sample of rows and columns (first image) or an uniform sample of pixels of the image (second one):

And this is what happens when you divide image into blocks and mix them randomly:

There is a funny and useful function called `jitter `which adds a small amount of noise to a numeric vector. What happens when you jitter every pixel of the image? As you can see, It becomes very vintage:

What if you transpose matrix? What if you change every color by another one? What if you change only a small range of them? What if you sum two images? What if you translate rgb colors into a grey scale? What if …? I answered some of these questions already and results are nice as well. After all, Marilyn can be represented as a simple matrix. Or maybe not.

```library("biOps")
library("abind")
#############################################################
#############################################################
plot(x)
#############################################################
#1. Swap
#############################################################
plot(imagedata(abind(x[,(ncol(x)/2):ncol(x),], x[,1:(ncol(x)/2),] , along=2)))
dev.copy(jpeg,filename="IMG01-Swap.jpg");
dev.off ();
#############################################################
#2. Artificial Symmetrical faces
#############################################################
plot(imagedata(abind(x[,1:(ncol(x)/2),], x[,(ncol(x)/2):1,] , along=2)))
dev.copy(jpeg,filename="IMG02-Symmetric1.jpg");
dev.off ();
plot(imagedata(abind(x[,ncol(x):(ncol(x)/2),], x[,(ncol(x)/2):ncol(x),] , along=2)))
dev.copy(jpeg,filename="IMG03-Symmetric2.jpg");
dev.off ();
#############################################################
#3. Uniform sampling over axis points
#############################################################
x2   <- aperm(array(255, dim = c(3, ncol(x), nrow(x))))
rows <- sample(1:nrow(x), round(nrow(x)*0.80), replace = FALSE)
cols <- sample(1:ncol(x), round(ncol(x)*0.80), replace = FALSE)
for (i in 1:length(rows))
{
for (j in 1: length(cols))
{
x2[rows[i], cols[j],1]<-x[rows[i], cols[j],1]
x2[rows[i], cols[j],2]<-x[rows[i], cols[j],2]
x2[rows[i], cols[j],3]<-x[rows[i], cols[j],3]
}
}
plot(imagedata(x2))
dev.copy(jpeg,filename="IMG04-Uniform1.jpg");
dev.off ();
#############################################################
#4. Uniform sampling over pixels
#############################################################
m2 <- matrix(rbinom(nrow(x)*ncol(x),1,0.5),nrow(x),ncol(x))
x4<- do.call(abind, c(list(x[,,1]*m2+(m2==0)*255,x[,,2]*m2+(m2==0)*255,x[,,3]*m2+(m2==0)*255), along = 3))
plot(imagedata(x4))
dev.copy(jpeg,filename="IMG05-Uniform2.jpg");
dev.off ();
#############################################################
#6. Jittering
#############################################################
x1<-mapply(as.matrix(x[,,1]), FUN=function(x)
{z<-round(x+jitter(0, amount=50))
if(z<0|z>255) x else z})
x1 <- matrix(x1, nrow = nrow(x),ncol = ncol(x))
x2<-mapply(as.matrix(x[,,2]), FUN=function(x)
{z<-round(x+jitter(0, amount=50))
if(z<0|z>255) x else z})
x2 <- matrix(x2, nrow = nrow(x),ncol = ncol(x))
x3<-mapply(as.matrix(x[,,3]), FUN=function(x)
{z<-round(x+jitter(0, amount=50))
if(z<0|z>255) x else z})
x3 <- matrix(x3, nrow = nrow(x),ncol = ncol(x))
x4<- do.call(abind, c(list(x1,x2,x3), along = 3))
plot(imagedata(x4))
dev.copy(jpeg,filename="IMG06-Jitter.jpg");
dev.off ();
#############################################################
#7. Mosaic
#############################################################
sptr <- 6 #Row splits
rnkr <- sample(1:sptr, size = sptr, replace = FALSE)
wthr <- floor(nrow(x)/sptr) #Splits width (row)
rnkr <- as.vector(sapply(rnkr, function (x) rep(x,wthr)))
rnkr <- rnkr*10E6+seq(1, length(rnkr), by=1)
rnkr <- rank(rnkr)
sptc <- round(ncol(x)/wthr)
rnkc <- sample(1:sptc, size = sptc, replace = FALSE)
wthc <- floor(ncol(x)/sptc) #Splits width (row)
rnkc <- as.vector(sapply(rnkc, function (x) rep(x,wthc)))
rnkc <- rnkc*10E6+seq(1, length(rnkc), by=1)
rnkc <- rank(rnkc)
x2<-x[1:length(rnkr),1:length(rnkc),]
x2<-x[rank(rnkr),rank(rnkc),]
plot(imagedata(x2))
dev.copy(jpeg,filename="IMG07-Mosaic.jpg");
dev.off ();
```