# Genetic Music: From Schoenberg to Bach

Bach, the epitome of a musician who strove all life long and finally acquired the ‘Habit of Perfection’, was a thoroughly imperfect human being (John Eliot Gardiner, Bach: Music in the Castle of Heaven)

Sometimes I dream awake and imagine I am a famous musician.  I fantasize being Paco de Lucía playing Mi niño Curro alone on the stage, Thom Yorke singing Fake plastic trees at Glastombury or Noel Gallagher singing Don’t look back in anger for a devoted crowd.

My parents gave me the opportunity to learn music, and this has been one of the best gifts I have received ever. I played the cello intensively until I had children but I still have enough skills to play some pieces. One of that is the Prelude of Suite No. 1 of J. S. Bach. It is very close to the limit of my possibilities but I love it. It is timeless, thrilling, provocative and elegant: an absolute masterpiece. I also imagine myself often playing it as well as my admired Yo-Yo Ma does.

The aim of this experiment is to discern first 4 beats of the prelude using a genetic algorithm. First of all, let’s listen our goal melody, created with tuneR package (sorry for the sound, Mr. Bach):

The frequency range of cello goes from 65.41 Hz to 987.77 Hz. Using the basic formula for the frequency of the notes, it means that a cello can produce 48 different notes. I generated the next codification for the 48 notes of the cello:

frequency (hz) note code
65.41 C2 a
69.30 C#2/Db2 b
73.42 D2 c
77.78 D#2/Eb2 d
82.41 E2 e
87.31 F2 f
92.50 F#2/Gb2 g
98.00 G2 h
103.83 G#2/Ab2 i
110.00 A2 j
116.54 A#2/Bb2 k
123.47 B2 l
130.81 C3 m
138.59 C#3/Db3 n
146.83 D3 o
155.56 D#3/Eb3 p
164.81 E3 q
174.61 F3 r
185.00 F#3/Gb3 s
196.00 G3 t
207.65 G#3/Ab3 u
220.00 A3 v
233.08 A#3/Bb3 w
246.94 B3 x
261.63 C4 y
277.18 C#4/Db4 z
293.66 D4 A
311.13 D#4/Eb4 B
329.63 E4 C
349.23 F4 D
369.99 F#4/Gb4 E
392.00 G4 F
415.30 G#4/Ab4 G
440.00 A4 H
466.16 A#4/Bb4 I
493.88 B4 J
523.25 C5 K
554.37 C#5/Db5 L
587.33 D5 M
622.25 D#5/Eb5 N
659.26 E5 O
698.46 F5 P
739.99 F#5/Gb5 Q
783.99 G5 R
830.61 G#5/Ab5 S
880.00 A5 T
932.33 A#5/Bb5 U
987.77 B5 V

So our goal melody is codified like this:

tAJHJAJAtAJHJAJAtCKJKCKCtCKJKCKCtEKJKEKEtEKJKEKEtFJHJFJFtFJHJFJF

I start with a population of 500 random melodies. All of them have 64 notes, the same length as the goal melody has. Given a melody, the algorithm compares it with the goal melody to calculate its fitness, with the following formula:

$fitness= {2}^{\displaystyle number of correct notes}$

For example, a melody with 5 correct notes has a fitness of 32. Being correct means being the right note in the right place. After measuring fitness of all melodies, I select 250 couples of individuals depending of its fitness (the more fitness, the more probability of being selected). Each couple generates two children for the next generation depending on certain probability, called crossover rate. Crossing operation is not always applied. Once two parents are selected, a random crossover point is chosen. At that point in both strings the genetic material from the left side of one parent is spliced to the material from the right side of other parent. The next figure illustrates the idea:

So two parents give birth to two children for the next generation. The last thing to do is mutate children. Once again, mutation is not always applied since it depends on a rate, usually small. Mutation introduces some new notes (new genetic material) to the next population. It increases convergence speed and reduces the probability to obtain a local optimum.

How many 32 -length melodies can be written with 48 notes? The answer is 4832, which is this extremely big number:

630.550.095.814.788.844.406.620.626.462.420.008.802.064.662.402.084.486

To understand how enormous is, let’s suppose we could work with Sunway TaihuLight, the fastest supercomputer in the world nowadays. This monster can do 93.000.000.000.000.000 floating-point operations per second so it will expend more than 214.995.831.974.513.789.322.026.202.008 years to calculate the fitness of all possible melodies: brute force is not an option.

A genetic algorithm does the job in just a few iterations. Best melodies introduce innovations which increase the average fitness of the whole population as well as its maximum fitness. Next table shows the evolution of an execution of the algorithm for a crossover rate equal of 75% and a mutation  rate of 1% (not exhaustive):

iteration best melody correct notes
1 OStxSTSbHwdsJAfTcRpoiNTRtRUxKhuRuKMcVNcBjRJNhENrVeFsPiegUpJHvRHw 7
5 tdbxSTSbHwdsJAfTcRpoiNTRtRITopoCPORzDdiFkEKrhEKtMHytiffzttJHvRHw 12
20 tAGHwdtUHzdMJATVACjJKVnetRQxKCKCtBKjqwiFkEKKhEKEMHyQiFfztUJHlRHF 25
50 tAJHwAJGjAJHJAJAtUCJKCkCtRUxKCKCtEKJKwKEtEKyhEKEMHyHrFfFtUJHJFHF 45
65 tAJHJAJGjAJHJAJAtUKJKCLCtCKxKCKCtEKJKwKEtEKyhEKEMHJHNFJFtFJHJFOF 52
80 tAJHJAJmtAJHJAJAtUKJKCLCtCKJKCKCtEKJKEKEtEKyMEKEMHJHJFJFtFJHJFOF 56
95 tAJHJAJjtAJHJAJAtUKJKCLCtCKJKCKCtEKJKEKEtEKJhEKEtFJHJFJFtFJHJFRF 59
110 tAJHJAJktAJHJAJAtUKJKCvCtCKJKCKCtEKJKEKEtEKJKEKEtFJHJFJFtFJHJFJF 61
125 tAJHJAJAtAJHJAJAtCKJKCKCtCKJKCKCtEKJKEKEtEKJKEKEtFJHJFJFtFJHJFJF 64

The optimum is reached in just 125 iterations. It is funny to merge the best melodies of some iterations. This sample blends four of them. The first one comes from the first initial population (the Schoenberg flavored) and the last one is our goal melody.  The other two were randomly picked from the rest iterations. It is nice to hear how the genetic algorithm turns randomness into the wonderful Bach’s melody:

This experiment was inspired by The Computational Beauty of Nature, a splendid book by Gary William Flake I strongly recommend you.

This is the code of the experiment:

library(tuneR)
library(stringdist)
library(dplyr)
#Function to calculate frequency
freq=function(n) 440*(2^(1/12))^n
#cello notes
notes=c("C2",
"C#2/Db2",
"D2",
"D#2/Eb2",
"E2",
"F2",
"F#2/Gb2",
"G2",
"G#2/Ab2",
"A2",
"A#2/Bb2",
"B2",
"C3",
"C#3/Db3",
"D3",
"D#3/Eb3",
"E3",
"F3",
"F#3/Gb3",
"G3",
"G#3/Ab3",
"A3",
"A#3/Bb3",
"B3",
"C4",
"C#4/Db4",
"D4",
"D#4/Eb4",
"E4",
"F4",
"F#4/Gb4",
"G4",
"G#4/Ab4",
"A4",
"A#4/Bb4",
"B4",
"C5",
"C#5/Db5",
"D5",
"D#5/Eb5",
"E5",
"F5",
"F#5/Gb5",
"G5",
"G#5/Ab5",
"A5",
"A#5/Bb5",
"B5")
#Table of frequencies
frequencies=data.frame(n=-33:14) %>%
mutate(frequency=round(freq(n),4),
note=notes,
code=c(letters, toupper(letters))[1:48])
#Codification of the goal melody
prelude="tAJHJAJAtAJHJAJAtCKJKCKCtCKJKCKCtEKJKEKEtEKJKEKEtFJHJFJFtFJHJFJF"
#Sample wav
if (exists("all_wave")) rm(all_wave)
frequencies %>%
filter(code==substr(prelude,1,1)) %>%
select(frequency) %>%
as.numeric %>%
sine(duration = 10000)->all_wave
for (i in 2:nchar(prelude))
frequencies %>%
filter(code==substr(prelude,i,i)) %>%
select(frequency) %>%
as.numeric %>%
sine(duration = 10000) %>% bind(all_wave, .)->all_wave
play(all_wave)
writeWave(all_wave, 'PreludeSample.wav')

popsize=500 #Population size
length=nchar(prelude)
genes=frequencies\$code
maxfitness=2^(1-(stringdist(prelude, prelude, method="hamming")-length))
maxiter=200 #Max number of iterations
iter=1
mutrate=0.01
#Initial population
replicate(popsize, sample(genes, length, replace = TRUE)) %>%
apply(2, function(x) paste(x,collapse="")) -> population
#Fitness evaluation
fitness=sapply(population, function(x) 2^(1-(stringdist(x, prelude, method="hamming")-length)), USE.NAMES=FALSE)
#Maximum fitness
maxfitenss_iter=max(fitness)
#Best melody
which((fitness)==max(fitness)) %>% min %>% population[.] ->bestfit
results=data.frame(iteration=iter, best_melody=bestfit, correct_notes=log(maxfitenss_iter, base = 2)-1)
#Execution of the algorithm
while(maxfitenss_iter<maxfitness & iter<maxiter)
{
population2=c()
for (i in 1:(popsize/2))
{
parents=sample(1:popsize, size=2, prob=fitness/sum(fitness), replace=FALSE)
mix=sample(1:(length-1), 1)

if (runif(1)>.25)
{
p1=paste0(substr(population[parents[1]],1,mix), substr(population[parents[2]],mix+1,length))
p2=paste0(substr(population[parents[2]],1,mix), substr(population[parents[1]],mix+1,length))
}
else
{
p1=population[parents[1]]
p2=population[parents[2]]
}
for (j in 1:length) if(runif(1)<mutrate) substr(p1,j,j)=sample(genes,1)
for (j in 1:length) if(runif(1)<mutrate) substr(p2,j,j)=sample(genes,1)
c(p1, p2) %>% c(population2)->population2
}
#New population
population=population2
fitness=sapply(population, function(x) 2^(1-(stringdist(x, prelude, method="hamming")-length)), USE.NAMES=FALSE)
which((fitness)==max(fitness)) %>% min %>% population[.] ->bestfit
print(paste0("Iteration ",iter, ": ", bestfit))
maxfitenss_iter=max(fitness)
iter=iter+1
data.frame(iteration=iter, best_melody=bestfit, correct_notes=log(maxfitenss_iter, base = 2)-1) %>% rbind(results) -> results
}

# The Sound Of Mandelbrot Set

Music is the pleasure the human soul experiences from counting without being aware that it is counting (Gottfried Leibniz)

I like the concept of sonification: translating data into sounds. There is a huge amount of contents in the Internet about this technique and there are several packages in R to help you to sonificate your data. Maybe one of the most accessible is tuneR, the one I choosed for this experiment. Do not forget to have a look to playitbyr: a package that allows you to listen to a data.frame in R by mapping columns onto sonic parameters, creating an auditory graph, as you can find in its website. It has a very similar syntaxis to ggplot. I will try to post something about playitbyr in the future.

Let me start plotting the Mandelbrot Set. I know you have seen it lot of times but it is very easy to plot in with R and result is extremely beautiful. Here you have four images corresponding to 12, 13, 14 and 15 iterations of the set’s generator. I like a lot how the dark blue halo around the Set evaporates as number of iterations increases.

And here you have the Set generated by 50 iterations. This is the main ingredient of the experiment:

Mandelbrot Set is generated by the recursive formula xt+1=xt2+c, with x0=0. A complex number c belongs to the Mandelbrot Set if its module after infinite iterations is finite. It is not possible to iterate a infinite number of times so every representation of Mandelbrot Set is just an approximation for a usually big amount of iterations. First image of Mandelbrot Set was generated in 1978 by Robert W. Brooks and Peter Matelski. You can find it here. I do not know how long it took to obtain it but you will spend only a couple of minutes to generate the ones you have seen before. It is amazing how computers have changed in this time!

This iterative equation is diabolical. To see just how pathological is, I transformed the succession of modules of xt generated by a given c in a succession of sounds. Since it is known that if one of this iterated complex numbers exceeds 2 in module then it is not in the Mandelbrot Set, frequencies of these sounds are bounded between 280 Hz (when module is equal to zero) and 1046 Hz (when module is equal or greater to 2). I called this function CreateSound. Besides the initial complex, you can choose how many notes and how long you want for your composition.

I tried with lot of numbers and results are funny. I want to stand out three examples from the rest:

• -1+0i gives the sequence 0, −1, 0, −1, 0 … which is bounded. Translated into music it sounds like an ambulance siren.
• -0.1528+1.0397i that is one of the generalized Feigenbaum points, around the Mandelbrot Set is conjetured to be self-similar. It sounds as a kind of Greek tonoi.
• -3/4+0.01i which presents a crazy slow divergence. I wrote a post some weeks ago about this special numbers around the neck of Mandelbrot Set and its relationship with PI.

All examples are ten seconds length. Take care with the size of the WAV file when you increase duration. You can create your own music files with the code below. If you want to download my example files, you can do it here. If you discover something interesting, please let me know.

Enjoy the music of Mandelbrot:

library(ggplot2)
library(reshape)
library(tuneR)
rm(list=ls())
# Create a grid of complex numbers
c.points <- outer(seq(-2.5, 1, by = 0.002),1i*seq(-1.5, 1.5, by = 0.002),'+')
z <- 0
for (k in 1:50) z <- z^2+c.points # Iterations of fractal's formula
c.points <- data.frame(melt(c.points))
colnames(c.points) <- c("r.id", "c.id", "point")
z.points <- data.frame(melt(z))
colnames(z.points) <- c("r.id", "c.id", "z.point")
mandelbrot <- merge(c.points, z.points, by=c("r.id","c.id")) # Mandelbrot Set
# Plotting only finite-module numbers
ggplot(mandelbrot[is.finite(-abs(mandelbrot\$z.point)), ], aes(Re(point), Im(point), fill=exp(-abs(z.point))))+
geom_tile()+theme(legend.position="none", axis.title.x = element_blank(), axis.title.y = element_blank())
#####################################################################################
# Function to translate numbers (complex modules) into sounds between 2 frequencies
#   the higher the module is, the lower the frequencie is
#   modules greater than 2 all have same frequencie equal to low.freq
#   module equal to 0 have high.freq
#####################################################################################
Module2Sound <- function (x, low.freq, high.freq)
{
if(x>2 | is.nan(x)) {low.freq} else {x*(low.freq-high.freq)/2+high.freq}
}
#####################################################################################
# Function to create wave. Parameters:
#    complex     : complex number to test
#    number.notes: number of notes to create (notes = iterations)
#    tot.duration.secs: Duration of the wave in seconds
#####################################################################################
CreateSound <- function(complex, number.notes, tot.duration.secs)
{
dur <- tot.duration.secs/number.notes
sep1 <- paste(", bit = 16, duration= ",dur, ", xunit = 'time'),sine(")
sep2 <- paste(", bit = 16, duration =",dur,",  xunit = 'time'))")
v.sounds <- c()
z <- 0
for (k in 1:number.notes)
{
z <- z^2+complex
v.sounds <- c(v.sounds, abs(z))
}
v.freqs <- as.vector(apply(data.frame(v.sounds), 1, FUN=Module2Sound, low.freq=280, high.freq=1046))
eval(parse(text=paste("bind(sine(", paste(as.vector(v.freqs), collapse = sep1), sep2)))
}
sound1 <- CreateSound(-3/4+0.01i     , 400 , 10) # Slow Divergence
sound2 <- CreateSound(-0.1528+1.0397i, 30  , 10) # Feigenbaum Point
sound3 <- CreateSound(-1+0i          , 20  , 10) # Ambulance Siren
writeWave(sound1, 'SlowDivergence.wav')
writeWave(sound2, 'FeigenbaumPoint.wav')
writeWave(sound3, 'AmbulanceSiren.wav')