**R – Fronkonstin**, and kindly contributed to R-bloggers)

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:

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 48^{32}, 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 |

35 | tAGHwAQUjAdsJAGAcUjJKCLCtRQxKCKCtEKAqwKEzEKJhEKEMHytIFfFtUJHJRHF | 35 |

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.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) % 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 }

**leave a comment**for the author, please follow the link and comment on their blog:

**R – Fronkonstin**.

R-bloggers.com offers

**daily e-mail updates**about R news and tutorials on topics such as: Data science, Big Data, R jobs, visualization (ggplot2, Boxplots, maps, animation), programming (RStudio, Sweave, LaTeX, SQL, Eclipse, git, hadoop, Web Scraping) statistics (regression, PCA, time series, trading) and more...