# Drawing 10 Million Points With ggplot: Clifford Attractors

For me, mathematics cultivates a perpetual state of wonder about the nature of mind, the limits of thoughts, and our place in this vast cosmos (Clifford A. Pickover – The Math Book: From Pythagoras to the 57th Dimension, 250 Milestones in the History of Mathematics)

I am a big fan of Clifford Pickover and I find inspiration in his books very often. Thanks to him, I discovered the harmonograph and the Parrondo’s paradox, among many other mathematical treasures. Apart of being a great teacher, he also invented a family of strange attractors wearing his name. Clifford attractors are defined by these equations:

$x_{n+1}\, =\, sin(a\, y_{n})\, +\, c\, cos(a\, x_{n}) \\ y_{n+1}\, =\, sin(b\, x_{n})\, +\, d\, cos(b\, y_{n}) \\$

There are infinite attractors, since a, b, c and d are parameters. Given four values (one for each parameter) and a starting point `(x0, y0)`, the previous equation defines the exact location of the point at step `n`, which is defined just by its location at `n-1`; an attractor can be thought as the trajectory described by a particle. This plot shows the evolution of a particle starting at `(x0, y0)=(0, 0)` with parameters `a=-1.24458046630025`, `b=-1.25191834103316`, `c=-1.81590817030519` and `d=-1.90866735205054` along 10 million of steps:

Changing parameters is really entertaining. Drawings have a sandy appearance:

From a technical point of view, the challenge is creating a data frame with all locations, since it must have 10 milion rows and must be populated sequentially. A very fast way to do it is using `Rcpp` package. To render the plot I use ggplot, which works quite well. Here you have the code to play with Clifford Attractors if you want:

```library(Rcpp)
library(ggplot2)
library(dplyr)

opt = theme(legend.position  = "none",
panel.background = element_rect(fill="white"),
axis.ticks       = element_blank(),
panel.grid       = element_blank(),
axis.title       = element_blank(),
axis.text        = element_blank())

cppFunction('DataFrame createTrajectory(int n, double x0, double y0,
double a, double b, double c, double d) {
// create the columns
NumericVector x(n);
NumericVector y(n);
x[0]=x0;
y[0]=y0;
for(int i = 1; i < n; ++i) {
x[i] = sin(a*y[i-1])+c*cos(a*x[i-1]);
y[i] = sin(b*x[i-1])+d*cos(b*y[i-1]);
}
// return a new data frame
return DataFrame::create(_["x"]= x, _["y"]= y);
}
')

a=-1.24458046630025
b=-1.25191834103316
c=-1.81590817030519
d=-1.90866735205054

df=createTrajectory(10000000, 0, 0, a, b, c, d)

png("Clifford.png", units="px", width=1600, height=1600, res=300)
ggplot(df, aes(x, y)) + geom_point(color="black", shape=46, alpha=.01) + opt
dev.off()
```

# Plants

Blue dragonflies dart to and fro
I tie my life to your balloon and let it go
(Warm Foothills, Alt-J)

In my last post I did some drawings based on L-Systems. These drawings are done sequentially. At any step, the state of the drawing can be described by the position (coordinates) and the orientation of the pencil. In that case I only used two kind of operators: drawing a straight line and turning a constant angle. Today I used two more symbols to do stack operations:

• “[“ Push the current state (position and orientation) of the pencil onto a pushdown
operations stack
• “]” Pop a state from the stack and make it the current state of the pencil (no line is drawn)

These operators allow to return to a previous state to continue drawing from there. Using them you can draw plants like these:

Each image corresponds to a different axiom, rules, angle and depth. I described these terms in my previous post. If you want to reproduce them you can find the code below (each image corresponds to a different set of axiom, rules, angle and depth parameters). Change colors, add noise to angles, try your own plants … I am sure you will find nice images:

```
library(gsubfn)
library(stringr)
library(dplyr)
library(ggplot2)

#Plant 1
axiom="F"
rules=list("F"="FF-[-F+F+F]+[+F-F-F]")
angle=22.5
depth=4

#Plant 2
axiom="X"
rules=list("X"="F[+X][-X]FX", "F"="FF")
angle=25.7
depth=7

#Plant 3
axiom="X"
rules=list("X"="F[+X]F[-X]+X", "F"="FF")
angle=20
depth=7

#Plant 4
axiom="X"
rules=list("X"="F-[[X]+X]+F[+FX]-X", "F"="FF")
angle=22.5
depth=5

#Plant 5
axiom="F"
rules=list("F"="F[+F]F[-F]F")
angle=25.7
depth=5

#Plant 6
axiom="F"
rules=list("F"="F[+F]F[-F][F]")
angle=20
depth=5

for (i in 1:depth) axiom=gsubfn(".", rules, axiom)

actions=str_extract_all(axiom, "\\d*\\+|\\d*\\-|F|L|R|\\[|\\]|\\|") %>% unlist

status=data.frame(x=numeric(0), y=numeric(0), alfa=numeric(0))
points=data.frame(x1 = 0, y1 = 0, x2 = NA, y2 = NA, alfa=90, depth=1)

for (action in actions)
{
if (action=="F")
{
x=points[1, "x1"]+cos(points[1, "alfa"]*(pi/180))
y=points[1, "y1"]+sin(points[1, "alfa"]*(pi/180))
points[1,"x2"]=x
points[1,"y2"]=y
data.frame(x1 = x, y1 = y, x2 = NA, y2 = NA,
alfa=points[1, "alfa"],
depth=points[1,"depth"]) %>% rbind(points)->points
}
if (action %in% c("+", "-")){
alfa=points[1, "alfa"]
points[1, "alfa"]=eval(parse(text=paste0("alfa",action, angle)))
}
if(action=="["){
data.frame(x=points[1, "x1"], y=points[1, "y1"], alfa=points[1, "alfa"]) %>%
rbind(status) -> status
points[1, "depth"]=points[1, "depth"]+1
}

if(action=="]"){
depth=points[1, "depth"]
points[-1,]->points
data.frame(x1=status[1, "x"], y1=status[1, "y"], x2=NA, y2=NA,
alfa=status[1, "alfa"],
depth=depth-1) %>%
rbind(points) -> points
status[-1,]->status
}
}

ggplot() +
geom_segment(aes(x = x1, y = y1, xend = x2, yend = y2),
lineend = "round",
colour="white",
data=na.omit(points)) +
coord_fixed(ratio = 1) +
theme(legend.position="none",
panel.background = element_rect(fill="black"),
panel.grid=element_blank(),
axis.ticks=element_blank(),
axis.title=element_blank(),
axis.text=element_blank())
```

# A Shiny App to Draw Curves Based on L-System

Don’t worry about a thing ’cause every little thing gonna be alright (Three Little Birds, Bob Marley)

One of my favourite books is The Computational Beauty of Nature by Gary William Flake where there is a fantastic chapter about fractals in which I discovered the L-Systems.

L-Systems were conceived  in 1968 by Aristide Lindenmayer, a Hungarian biologist, as a mathematical description of plant growth. Apart from the Wikipedia, there are many places on the Internet where you can read about them. If you are interested, don’t miss The Algorithmic Beauty of Plants, an awesome book by Przemysław Prusinkiewicz that you can obtain here for free.

Roughly speaking, a L-System is a very efficient way to make drawings. In its simplest way consists in two different actions: draw a straigh line and change the angle. This is just what you need, for example, to draw a square: draw a straigh line of  any length, turn 90 degrees (without drawing), draw another straigh line of the same length, turn 90 degrees in the same direction, draw, turn and draw again. Denoting `F` as the action of drawing a line of length d and `+` as turning 90 degrees right, the whole process to draw a square can be represented as `F+F+F+F`.

L-Systems are quite simple to program in R. You only need to substitute the rules iteratively into the axiom (I use `gsubfn` function to do it) and split the resulting chain into parts with `str_extract_all`, for example. The result is a set of very simple actions (draw or turn) that can be visualized with `ggplot` and its path geometry. There are four important parameters in L-Systems:

• The seed of the drawing, called axiom
• The substitutions to be applied iteratively, called rules
• How many times to apply substitutions, called depth
• Angle of each turning

For example, let’s define the next L-System:

• Axiom: `F-F-F-F`
• Rule: `F → F−F+F+FF−F−F+F`

The rule means that every `F` must be replaced by `F−F+F+FF−F−F+F` while `+` means right turning and `-` left one. After one iteration, the axiom is replaced by `F-F+F+FF-F-F+F-F-F+F+FF-F-F+F-F-F+F+FF-F-F+F-F-F+F+FF-F-F+F` and iterating again, the new string is `F-F+F+FF-F-F+F-F-F+F+FF-F-F+F+F-F+F+FF-F-F+F+F-F+F+FF-F-F+FF-F+F+FF-F-F+F-F-F+F+FF-F-F+F-F-F+F+FF-F-F+F+F-F+F+FF-F-F+F-F-F+F+FF-F-F+F-F-F+F+FF-F-F+F+F-F+F+FF-F-F+F+F-F+F+FF-F-F+FF-F+F+FF-F-F+F-F-F+F+FF-F-F+F-F-F+F+FF-F-F+F+F-F+F+FF-F-F+F-F-F+F+FF-F-F+F-F-F+F+FF-F-F+F+F-F+F+FF-F-F+F+F-F+F+FF-F-F+FF-F+F+FF-F-F+F-F-F+F+FF-F-F+F-F-F+F+FF-F-F+F+F-F+F+FF-F-F+F-F-F+F+FF-F-F+F-F-F+F+FF-F-F+F+F-F+F+FF-F-F+F+F-F+F+FF-F-F+FF-F+F+FF-F-F+F-F-F+F+FF-F-F+F-F-F+F+FF-F-F+F+F-F+F+FF-F-F+F`. As you can see, the length of the string grows exponentially. Converting last string into actions, produces this drawing, called Koch Island:

It is funny how different axioms and rules produce very different drawings. I have done a Shiny App to play with L-systems. Although it is quite simple, it has two interesting features I would like to undeline:

• Delay reactions with `eventReactive` to allow to set depth and angle values before refreshing the plot
• Build a dynamic UI that reacts to user input depending on the curve choosen

There are twelve curves in the application: Koch Island (and 6 variations), cuadratic snowflake, Sierpinsky triangle, hexagonal Gosper, quadratic Gosper and Dragon curve. These are their plots:

The definition of all these curves (axiom and rules) can be found in the first chapter of the Prusinkiewicz’s book. The magic comes when you modify angles and colors. These are some examples among the infinite number of possibilities that can be created:

I enjoyed a lot doing and playing with the app. You can try it here. If you do a nice drawing, please let me know in Twitter or dropping me an email. This is the code of the App:

`ui.R`:

```library(shiny)

shinyUI(fluidPage(
titlePanel("Curves based on L-systems"),

sidebarLayout(
sidebarPanel(
selectInput("cur", "Choose a curve:",
c("","Koch Island",
"Koch Variation 1",
"Koch Variation 2",
"Koch Variation 3",
"Koch Variation 4",
"Koch Variation 5",
"Koch Variation 6",
"Sierpinsky Triangle",
"Dragon Curve",
"Hexagonal Gosper Curve",
selected = ""),

conditionalPanel(
condition = "input.cur != \"\"",
uiOutput("Iterations")),

conditionalPanel(
condition = "input.cur != \"\"",
uiOutput("Angle")),

conditionalPanel(
condition = "input.cur != \"\"",
selectInput("lic", label = "Line color:", choices = colors(), selected = "black")),

conditionalPanel(
condition = "input.cur != \"\"",
selectInput("bac", label = "Background color:", choices = colors(), selected = "white")),

conditionalPanel(
condition = "input.cur != \"\"",
actionButton(inputId = "go", label = "Go!",
style="color: #fff; background-color: #337ab7; border-color: #2e6da4"))

),
mainPanel(plotOutput("curve", height="550px", width = "100%"))
)

))
```

`server.R`:

```library(shiny)
library(gsubfn)
library(stringr)
library(dplyr)
library(ggplot2)
library(rlist)

shinyServer(function(input, output) {

curves=list(
list(name="Koch Island",
axiom="F-F-F-F",
rules=list("F"="F-F+F+FF-F-F+F"),
angle=90,
n=2,
alfa0=90),
axiom="-F",
rules=list("F"="F+F-F-F+F"),
angle=90,
n=4,
alfa0=90),
list(name="Koch Variation 1",
axiom="F-F-F-F",
rules=list("F"="FF-F-F-F-F-F+F"),
angle=90,
n=3,
alfa0=90),
list(name="Koch Variation 2",
axiom="F-F-F-F",
rules=list("F"="FF-F-F-F-FF"),
angle=90,
n=4,
alfa0=90),
list(name="Koch Variation 3",
axiom="F-F-F-F",
rules=list("F"="FF-F+F-F-FF"),
angle=90,
n=3,
alfa0=90),
list(name="Koch Variation 4",
axiom="F-F-F-F",
rules=list("F"="FF-F--F-F"),
angle=90,
n=4,
alfa0=90),
list(name="Koch Variation 5",
axiom="F-F-F-F",
rules=list("F"="F-FF--F-F"),
angle=90,
n=5,
alfa0=90),
list(name="Koch Variation 6",
axiom="F-F-F-F",
rules=list("F"="F-F+F-F-F"),
angle=90,
n=4,
alfa0=90),
list(name="Sierpinsky Triangle",
axiom="R",
rules=list("L"="R+L+R", "R"="L-R-L"),
angle=60,
n=6,
alfa0=0),
list(name="Dragon Curve",
axiom="L",
rules=list("L"="L+R+", "R"="-L-R"),
angle=90,
n=10,
alfa0=90),
list(name="Hexagonal Gosper Curve",
axiom="L",
rules=list("L"="L+R++R-L--LL-R+", "R"="-L+RR++R+L--L-R"),
angle=60,
n=4,
alfa0=60),
axiom="-R",
rules=list("L"="LL-R-R+L+L-R-RL+R+LLR-L+R+LL+R-LR-R-L+L+RR-",
"R"="+LL-R-R+L+LR+L-RR-L-R+LRR-L-RL+L+R-R-L+L+RR"),
angle=90,
n=2,
alfa0=90))

output\$Iterations <- renderUI({ if (input\$cur!="") curve=list.filter(curves, name==input\$cur) else curve=list.filter(curves, name=="Koch Island") iterations=list.select(curve, n) %>% unlist
numericInput("ite", "Depth:", iterations, min = 1, max = (iterations+2))
})

output\$Angle <- renderUI({ curve=list.filter(curves, name==input\$cur) angle=list.select(curve, angle) %>% unlist
numericInput("ang", "Angle:", angle, min = 0, max = 360)
})

data <- eventReactive(input\$go, { curve=list.filter(curves, name==input\$cur) axiom=list.select(curve, axiom) %>% unlist
rules=list.select(curve, rules)[[1]]\$rules
alfa0=list.select(curve, alfa0) %>% unlist

for (i in 1:input\$ite) axiom=gsubfn(".", rules, axiom)
actions=str_extract_all(axiom, "\\d*\\+|\\d*\\-|F|L|R|\\[|\\]|\\|") %>% unlist

points=data.frame(x=0, y=0, alfa=alfa0)
for (i in 1:length(actions))
{
if (actions[i]=="F"|actions[i]=="L"|actions[i]=="R")
{
x=points[nrow(points), "x"]+cos(points[nrow(points), "alfa"]*(pi/180))
y=points[nrow(points), "y"]+sin(points[nrow(points), "alfa"]*(pi/180))
alfa=points[nrow(points), "alfa"]
points %>% rbind(data.frame(x=x, y=y, alfa=alfa)) -> points
}
else{
alfa=points[nrow(points), "alfa"]
points[nrow(points), "alfa"]=eval(parse(text=paste0("alfa",actions[i], input\$ang)))
}
}
return(points)
})

output\$curve <- renderPlot({
ggplot(data(), aes(x, y)) +
geom_path(color=input\$lic) +
coord_fixed(ratio = 1) +
theme(legend.position="none",
panel.background = element_rect(fill=input\$bac),
panel.grid=element_blank(),
axis.ticks=element_blank(),
axis.title=element_blank(),
axis.text=element_blank())
})

})
```

# Sunflowers for COLOURlovers

Andar, lo que es andar, anduve encima siempre de las nubes (Del tiempo perdido, Robe)

If you give importance to colours, maybe you know already COLOURlovers. As can be read in their website, COLOURlovers is a creative community where people from around the world create and share colors, palettes and patterns, discuss the latest trends and explore colorful articles… All in the spirit of love.

There is a R package called colourlovers which provides access to the COLOURlovers API. It makes very easy to choose nice colours for your graphics. I used clpalettes function to search for the top palettes of the website. Their names are pretty suggestive as well: Giant Goldfish, Thought Provoking, Adrift in Dreams, let them eat cake … Inspired by this post I have done a Shiny app to create colored flowers using that palettes. Seeds are arranged according to the golden angle. One example:

Some others:

You can play with the app here.

If you want to do your own sunflowers, here you have the code. This is the `ui.R` file:

```library(colourlovers)
library(rlist)
top=clpalettes('top')
sapply(1:length(top), function(x) list.extract(top, x)\$title)-&gt;titles

fluidPage(
titlePanel("Sunflowers for COLOURlovers"),
fluidRow(
column(3,
wellPanel(
selectInput("pal", label = "Palette:", choices = titles),
sliderInput("nob", label = "Number of points:", min = 200, max = 500, value = 400, step = 50)
)
),
mainPanel(
plotOutput("Flower")
)
)
)
```

And this is the `server.R` one:

```library(shiny)
library(ggplot2)
library(colourlovers)
library(rlist)
library(dplyr)

top=clpalettes('top')
sapply(1:length(top), function(x) list.extract(top, x)\$title)->titles

CreatePlot = function (ang=pi*(3-sqrt(5)), nob=150, siz=15, sha=21, pal="LoversInJapan") {

list.extract(top, which(titles==pal))\$colors %>%
unlist %>%
as.vector() %>%
paste0("#", .) -> all_colors

colors=data.frame(hex=all_colors, darkness=colSums(col2rgb(all_colors)))
colors %>% arrange(-darkness)->colors

background=colors[1,"hex"] %>% as.character

colors %>% filter(hex!=background) %>% .[,1] %>% as.vector()->colors

ggplot(data.frame(r=sqrt(1:nob), t=(1:nob)*ang*pi/180), aes(x=r*cos(t), y=r*sin(t)))+
geom_point(colour=sample(colors, nob, replace=TRUE, prob=exp(1:length(colors))), aes(size=(nob-r)), shape=16)+
scale_x_continuous(expand=c(0,0), limits=c(-sqrt(nob)*1.4, sqrt(nob)*1.4))+
scale_y_continuous(expand=c(0,0), limits=c(-sqrt(nob)*1.4, sqrt(nob)*1.4))+
theme(legend.position="none",
panel.background = element_rect(fill=background),
panel.grid=element_blank(),
axis.ticks=element_blank(),
axis.title=element_blank(),
axis.text=element_blank())}

function(input, output) {
output\$Flower=renderPlot({
CreatePlot(ang=180*(3-sqrt(5)), nob=input\$nob, siz=input\$siz, sha=as.numeric(input\$sha), pal=input\$pal)
}, height = 550, width = 550 )}
```

# The Ex Libris Generator

Go ahead stomp your feet on the floorboards
Clap your hands if that’s really what you came here for
(Heaven, The Milk Carton Kids)

Inspired by curves created by the harmonograph, I have done a Shiny app to generate random images that you can personalize and use as an Exlibris.  You can try the App here. For me, an exlibris (also known as bookplates) can be a nice, original and useful present for book-lovers. This is an example:

More examples:

I always put the code at the end of my posts. Since I always have doubts about how many people are interested in what I do, today will be different. I will share the code with those who ask it to me in any of the following ways:

• Sending me a direct message on Twitter
• Droping me an email

Cheers!

# Three Shiny Apps to Celebrate the Beauty of Maths

Mathematics knows no races or geographic boundaries; for mathematics, the cultural world is one country (David Hilbert)

One of the best decisions I took this year related with this blog was to move it to my own self-hosted domain using WordPress.org. It allows to me, for example, to embed dynamic JavaScript visualizations like this one. Another thing I can do now is to upload my Shiny Apps to share them with my readers. In this post I have gathered three Apps I made some time ago; you can play with them as well as get the code I wrote for each one:

• The Harmonograph: This App simulates harmonograph drawings. An harmonograph is a mechanism which draws trajectories by means of two pendulums: one moves a pencil and the other one moves a platform with a piece of paper on it. Click here to try it.
• Shiny Wool Skeins: This App, inspired by this post, creates a plot consisting of chords inside a circle . You can change colors as well as the number and quality of the chords. Click here to try it.
• The Coaster Maker: With this App you can create your own coasters using hypocicloids. Click here to try it.

I want to thank to my friend Jorge, without whom I would not have been able to make Shiny work in my server.

# Amazing Things That Happen When You Toss a Coin 12 Times

If there is a God, he’s a great mathematician (Paul Dirac)

Imagine you toss a coin 12 times and you count how many heads and tails you are obtaining after each throwing (the coin is equilibrated so the probability of head or tail is the same). At some point, it can happen that number of heads and number of tails are the same. For example, if you obtain the sequence T-H-T-T-H-T-H-H-T-T-H-H, after the second throwing, number of heads is equal to number of tails (and both equal to one). It happens again after the 8th throwing and after last one. In this example, the last throwing where equallity occurs is the number 12. Obviously, equallity can only be observed in even throwings.

If you repeat the experiment 10.000 times you will find something like this if you draw the relative frequency of the last throwing where cumulated number of heads is equal to the one of tails:

From my point of view there are three amazing things in this plot:

1. It is symmetrical, so `prob(n)=prob(12-n)`
2. The least likely throwing to obtain the last equality is the central one.
3. As a corollary, the most likely is not obtaining any equality (number of heads never are the same than number of tails) or obtaining last equality in the last throwing: two extremely different scenarios with the same chances to be observed.

Behind the simplicity of tossing coins there is a beautiful universe of mathematical surprises.

```library(dplyr)
library(ggplot2)
library(scales)
tosses=12
iter=10000
results=data.frame(nmax=numeric(0), count=numeric(0), iter=numeric(0))
tmp=data.frame(nmax=numeric(0))
for (j in 1:iter)
{
data.frame(x=sample(c(-1,1), size=tosses, replace=TRUE)) %>%
add_rownames(var = "n") %>%
mutate(cumsum = cumsum(x)) %>% filter(cumsum==0) %>%
summarize(nmax=max(as.numeric(n))) %>% rbind(tmp)->tmp
}
tmp %>%
group_by(nmax) %>%
summarize(count=n()) %>%
mutate(nmax=ifelse(is.finite(nmax), nmax, 0), iter=iter) %>%
rbind(results)->results
opts=theme(
panel.background = element_rect(fill="darkolivegreen1"),
panel.border = element_rect(colour="black", fill=NA),
axis.line = element_line(size = 0.5, colour = "black"),
axis.ticks = element_line(colour="black"),
panel.grid.major = element_line(colour="white", linetype = 1),
panel.grid.minor = element_blank(),
axis.text.y = element_text(colour="black"),
axis.text.x = element_text(colour="black"),
text = element_text(size=20),
legend.key = element_blank(),
plot.title = element_text(size = 30)
)
ggplot(results, aes(x=nmax, y=count/iter)) +
geom_line(size=2, color="green4")+
geom_point(size=8, fill="green4", colour="darkolivegreen1",pch=21)+
scale_x_continuous(breaks = seq(0, tosses, by=2))+
scale_y_continuous(labels=percent, limits=c(0, .25))+
labs(title="What happens when you toss a coin 12 times?",
x="Last throwing where cumulated #tails = #heads",
y="Probability (estimated)")+opts
```

# The Unbereable Insolence of Prime Numbers or (Playing to be Ulam)

So rock me mama like a wagon wheel, rock me mama anyway you feel (Wagon Wheel, Old Crow Medicine Show)

This is the third iteration of Hilbert curve. I placed points in its corners. Since the curve has beginning and ending, I labeled each vertex with the order it occupies:Dark green vertex are those labeled with prime numbers and light ones with non-prime. This is the sixth iteration colored as I described before (I removed lines and labels):

Previous plot has 4.096 points. There are 564 primes lower than 4.096. What If I color 564 points randomly instead coloring primes? This is an example:

Do you see any difference? I do. Let me place both images together (on the left, the one with primes colored):

The dark points are much more ordered in the first plot. The second one is more noisy. This is my particular tribute to Stanislaw Ulam and its spiral: one of the most amazing fruits of boredom in the history of mathematics.

This is the code:

```library(reshape2)
library(dplyr)
library(ggplot2)
library(pracma)
opt=theme(legend.position="none",
panel.background = element_rect(fill="white"),
panel.grid=element_blank(),
axis.ticks=element_blank(),
axis.title=element_blank(),
axis.text=element_blank())
hilbert = function(m,n,r) {
for (i in 1:n)
{
tmp=cbind(t(m), m+nrow(m)^2)
m=rbind(tmp, (2*nrow(m))^r-tmp[nrow(m):1,]+1)
}
melt(m) %>% plyr::rename(c("Var1" = "x", "Var2" = "y", "value"="order")) %>% arrange(order)}
iter=3 #Number of iterations
df=hilbert(m=matrix(1), n=iter, r=2)
subprimes=primes(nrow(df))
df %>%  mutate(prime=order %in% subprimes,
random=sample(x=c(TRUE, FALSE), size=nrow(df), prob=c(length(subprimes),(nrow(df)-length(subprimes))), replace = TRUE)) -> df
#Labeled (primes colored)
ggplot(df, aes(x, y, colour=prime)) +
geom_path(color="gray75", size=3)+
geom_point(size=28)+
scale_colour_manual(values = c("olivedrab1", "olivedrab"))+
scale_x_continuous(expand=c(0,0), limits=c(0,2^iter+1))+
scale_y_continuous(expand=c(0,0), limits=c(0,2^iter+1))+
geom_text(aes(label=order), size=8, color="white")+
opt
#Non labeled (primes colored)
ggplot(df, aes(x, y, colour=prime)) +
geom_point(size=5)+
scale_colour_manual(values = c("olivedrab1", "olivedrab"))+
scale_x_continuous(expand=c(0,0), limits=c(0,2^iter+1))+
scale_y_continuous(expand=c(0,0), limits=c(0,2^iter+1))+
opt
#Non labeled (random colored)
ggplot(df, aes(x, y, colour=random)) +
geom_point(size=5)+
scale_colour_manual(values = c("olivedrab1", "olivedrab"))+
scale_x_continuous(expand=c(0,0), limits=c(0,2^iter+1))+
scale_y_continuous(expand=c(0,0), limits=c(0,2^iter+1))+
opt
```

# Going Bananas With Hilbert

It seemed that everything is in ruins, and that all the basic mathematical concepts have lost their meaning (Naum Vilenkin, Russian mathematician, regarding to the discovery of Peano’s curve)

Giuseppe Peano found in 1890 a way to draw a curve in the plane that filled the entire space: just a simple line covering completely a two dimensional plane. Its discovery meant a big earthquake in the traditional structure of mathematics. Peano’s curve was the first but not the last: one of these space-filling curves was discovered by Hilbert and takes his name. It is really beautiful:

Hilbert’s curve can be created iteratively. These are the first six iterations of its construction:

As you will see below, R code to create Hilbert’s curve is extremely easy. It is also very easy to play with the curve, altering the order in which points are sorted. Changing the initial `matrix(1)` by some other number, resulting curves are quite appealing:

Let’s go futher. Changing `ggplot` geometry from `geom_path` to `geom_polygon` generate some crazy pseudo-tessellations:

And what if you change the matrix exponent?

And what if you apply polar coordinates?

We started with a simple line and with some small changes we have created fantastical images. And all these things only using black and white. Do you want to add some colors? Try with the following code (if you draw something interesting, please let me know):

```library(reshape2)
library(dplyr)
library(ggplot2)
opt=theme(legend.position="none",
panel.background = element_rect(fill="white"),
panel.grid=element_blank(),
axis.ticks=element_blank(),
axis.title=element_blank(),
axis.text=element_blank())
hilbert = function(m,n,r) {
for (i in 1:n)
{
tmp=cbind(t(m), m+nrow(m)^2)
m=rbind(tmp, (2*nrow(m))^r-tmp[nrow(m):1,]+1)
}
melt(m) %>% plyr::rename(c("Var1" = "x", "Var2" = "y", "value"="order")) %>% arrange(order)}
# Original
ggplot(hilbert(m=matrix(1), n=1, r=2), aes(x, y)) + geom_path()+ opt
ggplot(hilbert(m=matrix(1), n=2, r=2), aes(x, y)) + geom_path()+ opt
ggplot(hilbert(m=matrix(1), n=3, r=2), aes(x, y)) + geom_path()+ opt
ggplot(hilbert(m=matrix(1), n=4, r=2), aes(x, y)) + geom_path()+ opt
ggplot(hilbert(m=matrix(1), n=5, r=2), aes(x, y)) + geom_path()+ opt
ggplot(hilbert(m=matrix(1), n=6, r=2), aes(x, y)) + geom_path()+ opt
# Changing order
ggplot(hilbert(m=matrix(.5), n=5, r=2), aes(x, y)) + geom_path()+ opt
ggplot(hilbert(m=matrix(0), n=5, r=2), aes(x, y)) + geom_path()+ opt
ggplot(hilbert(m=matrix(tan(1)), n=5, r=2), aes(x, y)) + geom_path()+ opt
ggplot(hilbert(m=matrix(3), n=5, r=2), aes(x, y)) + geom_path()+ opt
ggplot(hilbert(m=matrix(-1), n=5, r=2), aes(x, y)) + geom_path()+ opt
ggplot(hilbert(m=matrix(log(.1)), n=5, r=2), aes(x, y)) + geom_path()+ opt
ggplot(hilbert(m=matrix(-15), n=5, r=2), aes(x, y)) + geom_path()+ opt
ggplot(hilbert(m=matrix(-0.001), n=5, r=2), aes(x, y)) + geom_path()+ opt
# Polygons
ggplot(hilbert(m=matrix(log(1)), n=4, r=2), aes(x, y)) + geom_polygon()+ opt
ggplot(hilbert(m=matrix(.5), n=4, r=2), aes(x, y)) + geom_polygon()+ opt
ggplot(hilbert(m=matrix(tan(1)), n=5, r=2), aes(x, y)) + geom_polygon()+ opt
ggplot(hilbert(m=matrix(-15), n=4, r=2), aes(x, y)) + geom_polygon()+ opt
ggplot(hilbert(m=matrix(-25), n=4, r=2), aes(x, y)) + geom_polygon()+ opt
ggplot(hilbert(m=matrix(0), n=4, r=2), aes(x, y)) + geom_polygon()+ opt
ggplot(hilbert(m=matrix(1000000), n=4, r=2), aes(x, y)) + geom_polygon()+ opt
ggplot(hilbert(m=matrix(-1), n=4, r=2), aes(x, y)) + geom_polygon()+ opt
ggplot(hilbert(m=matrix(-.00001), n=4, r=2), aes(x, y)) + geom_polygon()+ opt
# Changing exponent
gplot(hilbert(m=matrix(log(1)), n=4, r=-1), aes(x, y)) + geom_polygon()+ opt
ggplot(hilbert(m=matrix(.5), n=4, r=-2), aes(x, y)) + geom_polygon()+ opt
ggplot(hilbert(m=matrix(tan(1)), n=4, r=6), aes(x, y)) + geom_polygon()+ opt
ggplot(hilbert(m=matrix(-15), n=3, r=sin(2)), aes(x, y)) + geom_polygon()+ opt
ggplot(hilbert(m=matrix(-25), n=4, r=-.0001), aes(x, y)) + geom_polygon()+ opt
ggplot(hilbert(m=matrix(0), n=4, r=200), aes(x, y)) + geom_polygon()+ opt
ggplot(hilbert(m=matrix(1000000), n=3, r=.5), aes(x, y)) + geom_polygon()+ opt
ggplot(hilbert(m=matrix(-1), n=4, r=sqrt(2)), aes(x, y)) + geom_polygon()+ opt
ggplot(hilbert(m=matrix(-.00001), n=4, r=52), aes(x, y)) + geom_polygon()+ opt
# Polar coordinates
ggplot(hilbert(m=matrix(1), n=4, r=2), aes(x, y)) + geom_polygon()+ coord_polar()+opt
ggplot(hilbert(m=matrix(-1), n=5, r=2), aes(x, y)) + geom_polygon()+ coord_polar()+opt
ggplot(hilbert(m=matrix(.1), n=2, r=.5), aes(x, y)) + geom_polygon()+ coord_polar()+opt
ggplot(hilbert(m=matrix(1000000), n=2, r=.1), aes(x, y)) + geom_polygon()+ coord_polar()+opt
ggplot(hilbert(m=matrix(.25), n=3, r=3), aes(x, y)) + geom_polygon()+ coord_polar()+opt
ggplot(hilbert(m=matrix(tan(1)), n=5, r=1), aes(x, y)) + geom_polygon()+ coord_polar()+opt
ggplot(hilbert(m=matrix(1), n=4, r=1), aes(x, y)) + geom_polygon()+ coord_polar()+opt
ggplot(hilbert(m=matrix(log(1)), n=3, r=sin(2)), aes(x, y)) + geom_polygon()+ coord_polar()+opt
ggplot(hilbert(m=matrix(-.0001), n=4, r=25), aes(x, y)) + geom_polygon()+ coord_polar()+opt
```

# Polar Circles

You cannot find peace by avoiding life (Virginia Woolf)

Combining polar coordinates, `RColorBrewer` palettes, `ggplot2` and a simple trigonometric function to define the width of the tiles is easy to produce nice circular plots like these:

Do you want to try? Here you have the code:

```library(ggplot2)
library(dplyr)
library(RColorBrewer)
n=500
m=50
w=sapply(seq(from=-3.5*pi, to=3.5*pi, length.out=n), function(x) {abs(sin(x))})
x=c(1)
for (i in 2:n) {x[i]=x[i-1]+1/2*(w[i-1]+w[i])}
expand.grid(x=x, y=1:m) %>%
mutate(w=rep(w, m))-> df
opt=theme(legend.position="none",
panel.background = element_rect(fill="white"),
panel.grid=element_blank(),
axis.ticks=element_blank(),
axis.title=element_blank(),
axis.text=element_blank())
ggplot(df, aes(x=x,y=y))+geom_tile(aes(fill=x, width=w))+
scale_fill_gradient(low=brewer.pal(9, "Greens")[1], high=brewer.pal(9, "Greens")[9])+
coord_polar(start = runif(1, min = 0, max = 2*pi))+opt
ggplot(df, aes(x=x,y=y))+geom_tile(aes(fill=w, width=w))+
scale_fill_gradient(low=brewer.pal(9, "Reds")[1], high=brewer.pal(9, "Reds")[9])+
coord_polar(start = runif(1, min = 0, max = 2*pi))+opt
ggplot(df, aes(x=x,y=y))+geom_tile(aes(fill=y, width=w))+
scale_fill_gradient(low=brewer.pal(9, "Purples")[1], high=brewer.pal(9, "Purples")[9])+
coord_polar(start = runif(1, min = 0, max = 2*pi))+opt
ggplot(df, aes(x=x,y=y))+geom_tile(aes(fill=w*y, width=w))+
scale_fill_gradient(low=brewer.pal(9, "Blues")[9], high=brewer.pal(9, "Blues")[1])+
coord_polar(start = runif(1, min = 0, max = 2*pi))+opt
```