index.Rmd

---
title: "Markovian Blues"
subtitle: "Analysis and Generation of 12-Bar Blues Melodies Using VLMCs"
author: "Lance J. Fernando"
date: "4/20/2018"
output: 
  html_document:
    toc: true # table of content true
    number_sections: true  ## if you want number sections at each table header
    theme: spacelab  # many options for theme, this one is my favorite.
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo=TRUE, messages=FALSE, warnings=FALSE, eval=FALSE)
setwd("~/Desktop/USFSpring2018/Stochastic_Processes/project/")
```

```{r load}
library(PST)
library(dplyr)
library(ggplot2)
library(magrittr)
library(reticulate)
library(audio)
```

# Song Preprocessing
```{r}
# Transposes a vector of pitches to the key of
# C. Accidentals will be denoted with flats 
# instead of sharp.
#
# args:
#     pitches - Vector of melodic pitches
#     key - Original key of given melodic pitches
#
# returns:
#     transposed vector with flatted accidentals
#
transposeC <- function(pitches, key){
  all_pitches <- c("C", "C#", "Db", "D", "D#",
                   "Eb", "E", "F", "F#", "Gb",
                   "G", "G#", "Ab", "A", "A#",
                   "Bb", "B")
  flat_pitches <- c("C", "Db", "D", "Eb", "E", "F",
                    "Gb", "G", "Ab", "A", "Bb", "B")
  
  # If the pitch is notated as sharp, return the flatted version
  get_flat <- function(pitch){
    if(substr(pitch, 2, 2) == "#"){
        return(all_pitches[which(all_pitches == pitch) + 1])
    }else{
      return(pitch)
    }
  }
  
  # Flatted version of pitches
  flat <- sapply(X = pitches, FUN = get_flat) %>% unname()
  
  # Flatted key (if initially indicated as sharp)
  flat_key <- get_flat(key)
  
  key_idx <- which(flat_pitches == flat_key)
  
  # Calculating the number of half steps to move to C
  note_diff <- key_idx - 1
  
  # Computing the transposed pitch vector
  transposed <- sapply(X = flat, FUN = function(p){
    if(p == "R"){
      return(p)
    }else{
      
      p_idx <- which(flat_pitches == p)
      if(p_idx <= note_diff){
        new_p <- flat_pitches[(12 + p_idx) - note_diff]
        return(new_p)
        
      }else{
        new_p <- flat_pitches[p_idx - note_diff]
        return(new_p)
      }
        
    }
  }) %>% unlist() %>% unname()
  
  return(transposed)
}

# Creates a column of sections given a vector of 
# note duration values. The total sum of durations
# must be 12 for a 12-bar-blues progression. An error
# is thrown if this sum is violated. 
#
# args:
#     durations - vector of note values for each note in a melody
#     structure - the type of chord structure to follow 
#                 (default: 12-bar-blues which groups bars in sets of 2s)
#                   
getSections <- function(durations, structure = c("Ia", "Ib", "IV", "Ic", "V", "Id")){
  if(sum(durations) != 12){
    
    return(cat("ERROR \n Duration sum: ", sum(durations), "\n",
               "Duration sum must be 12"))
  }
  sections <- character(0)
  cur_sec_idx <- 1
  cur_sec <- structure[cur_sec_idx]
  cur_sum <- 0
  
  # Threshold to classify next section
  sec_size <- 12/length(structure)
  
  for(i in 1:length(durations)){
    cur_sum <- cur_sum + durations[i]
    sections <- c(sections, cur_sec)
    
    # Can't calculate exact 0, so we check if difference
    # is smaller than a sixteenth note (which is the smallest note dur used)
    if(abs(cur_sum - sec_size) < 0.0625){
      cur_sec_idx <- cur_sec_idx + 1
      cur_sec <- structure[cur_sec_idx]
      cur_sum <- 0
    }else{
      next
    }
  }
  return(sections)
}


# Creates a column of sections given a vector of 
# note duration values. The total sum of durations
# must be 12 for a 12-bar-blues progression. An error
# is thrown if this sum is violated. 
#
# Takes a set of vectors specifying notes of a song
# and writes it to a csv file after transposing
# to C and converting note durations to categorical
# note symbols.
#
# args:
#     title - title of song used to save to csv
#     key - current key of song
#     n_pitches - vector of all the pitches in the melody
#     n_durr - numeric vector indicating the duration of each pitch in n_pitches
#
# output:
#     CSV file (pitch string, duration float, value string, section string)
writeSong <- function(title, key, n_pitches, n_durr){
  transpose <- transposeC(pitches = n_pitches, key = key)
  sections <- getSections(n_durr)
  
  value <- rep(NA, length(n_durr))
  value[n_durr == 1/16] <- "s"
  value[n_durr == 1/12] <- 't'
  value[n_durr == 1/8] <- 'e'
  value[n_durr == 1/4] <- 'q'
  value[n_durr == 3/8] <- 'dq'
  value[n_durr == 2/4] <- 'h'
  value[n_durr == 5/8] <- 'he'
  value[n_durr == 3/4] <- 'dh'
  value[n_durr == 7/8] <- 'ddh'
  value[n_durr == 1] <- 'w'
  
  song <- data.frame(pitch = transpose,
                     duration = n_durr,
                     value = value,
                     section = sections)

  write.csv(song, paste("songs/", title, ".csv", sep=""),
            row.names = FALSE)
}
```

```{r convert}
# Takes a given song in the form of a dataframe
# and breaks the notes into 1/8 time-value states.
# This results in all twelve-bar blues songs having
# 8*12 = 96 observations as opposed to one observation
# for each note. The new note values will indicate the
# pitch as well as an _S or _L symbol indicating whether 
# the current note stays or leaves. A Bb quarter-note will 
# have observations c(Bb_S, Bb_L) having duration 1/8 + 1/8 = 1/4
# whereas a Bb dotted-half will have c(Bb_S, Bb_S, Bb_S, Bb_S, Bb_S, Bb_L)
# with a total duration of 1/8 * 6 = 3/4. 
#
# args:
#     df - Given dataframe created by writeSong() function
#
# returns:
#     new dataframe (w/ dim:96 x 2) with two columns :
#         (state string, section string)
recodeSongs <- function(df){
  new_vals <- character(0)
  triplet_num <- 1
  for(i in 1:nrow(df)){
    value <- df$value[i]
    pitch <- df$pitch[i]
    
    
    if(value == 'e'){
      new_vals <- c(new_vals, paste(pitch,"L", sep="_"))
    }else if(value == "q"){
      new_vals <- c(new_vals, paste(pitch,"S", sep="_"),paste(pitch,"L",sep="_"))
    }else if(value == "dq"){
      new_vals <- c(new_vals, rep(paste(pitch,"S", sep="_"), 2),
                    paste(pitch,"L", sep="_"))
    }else if(value == "h"){
      new_vals <- c(new_vals, rep(paste(pitch,"S", sep="_"), 3),
                    paste(pitch,"L",sep="_"))
    }else if(value == "he"){
      new_vals <- c(new_vals, rep(paste(pitch,"S", sep="_"), 4),
                    paste(pitch,"L", sep="_"))
    }else if(value == "dh"){
      new_vals <- c(new_vals, rep(paste(pitch,"S", sep="_"), 5),
                    paste(pitch,"L", sep="_"))
    }else if(value == "ddh"){
      new_vals <- c(new_vals, rep(paste(pitch,"S", sep="_"), 6),
                    paste(pitch,"L", sep="_"))
    }else if(value == "w"){
      new_vals <- c(new_vals, rep(paste(pitch,"S", sep="_"), 7),
                    paste(pitch,"L", sep="_"))
    }else if(value == "t"){

      if(triplet_num == 1){
        pitch2 <- df$pitch[i + 1]
        new_vals <- c(new_vals, 
                      paste(pitch,"S", sep="_"),paste(pitch2, "L", sep="_"))
        triplet_num <- triplet_num + 1
        
      }else if(triplet_num == 2){
        triplet_num <- triplet_num + 1
      }else{
        triplet_num <- 1
      }
    }else{
      print("Unknown note duration")
      break
    }
 
  }
  new_secs <- c(rep("Ia", 16), rep("Ib", 16), rep("IV", 16), rep("Ic", 16), rep("V", 16), rep("Id",16))
  return(data.frame(state = new_vals,
                    section = new_secs,
                    stringsAsFactors = FALSE))
}

# Converts a granular version of notes broken down
# by 1/8 segments (i.e., Bb_S - Bb_L) into regular 
# pitches and values (i.e., pitch:Bb, value:q).
# 
# args:
#     genSeq - single vector output from generate() containing 1/8 segment breakdown
#
# returns:
#     dataframe(pitch, value)
convertSeqToNotes <- function(genSeq){
  
  new_seq <- unlist(genSeq) %>% as.character()

  # Replacing last character with terminating pitch
  term_pitch <- sub(pattern = "_S", replacement = "_L", x = new_seq[length(new_seq)])
  new_seq[length(new_seq)] <- term_pitch
  
  
  new_pitch <- character(0)
  new_val <- character(0)
  midi_val <- numeric(0)
  #Duration accumulator
  dur <- 0

  for(i in 1:length(new_seq)){
      
    symbol <- new_seq[i]
    if(regexpr(pattern = "L", text = symbol) == -1){ # If it's a held note
      dur <- dur + 1
    }else{# If it's an end note
      dur <- dur + 1
      new_pitch <- c(new_pitch, sub(pattern = "_L", replacement ="", x = symbol))
      new_val <- c(new_val, switch(dur,
                                   `1` = 'e',
                                   `2` = 'q',
                                   `3` = 'dq',
                                   `4` = 'h',
                                   `5` = 'he',
                                   `6` = 'dh',
                                   `7` = 'ddh',
                                   `8` = 'w'))
      midi_val <- c(midi_val, dur)
      dur <- 0
    }
  }
  
  midi_pitch <- recode(new_pitch,  C=60, Db=61, D=62, Eb=63, E=64, F=65,
         Gb=66, G=67, Ab=68, A=69, A=70, Bb=71, B=72, R=0) %>% as.numeric()

  return(data.frame(pitch = new_pitch,
                    mpitch = midi_pitch,
                    value = new_val,
                    mval = midi_val,
                    stringsAsFactors = FALSE))
}
```


# Example Writing Songs To CSV
Below is an example of writing *Creole Love Song* into CSV format. Notating this is obviously a tedious process and further examination into automating this would benefit this project greatly! This chunk of code will not be run as all the data was already created and stored in the */songs/* directory.
```{r writeSongs, eval=FALSE}
title <- "Creole_Love_Song"
key <- "C"
n_pitch <- c('G',
             'R', 'E', 'G', 'E', 'G', 'E', 'G', 'E',
             'G',
             'R', 'E', 'G', 'E', 'G', 'E', 'G', 'E',
             'A',
             'R', 'F', 'A', 'F', 'A', 'F', 'A', 'F',
             'A', 'A', 'G',
             'G', 'R', 'G', 'A', 'A#',
             'B', 'A', 'G', 'F',
             'C', 'A', 'G', 'A','G', 'Eb',
             'C',
             'R', 'E', 'G', 'E', 'G', 'E', 'G', 'E')
n_durr <- c(1,
            1/8, 1/8, 1/8, 1/8, 1/8, 1/8, 1/8, 1/8,
            1,
            1/8, 1/8, 1/8, 1/8, 1/8, 1/8, 1/8, 1/8,
            1,
            1/8, 1/8, 1/8, 1/8, 1/8, 1/8, 1/8, 1/8,
            3/4, 1/8, 1/8,
            1/2, 1/4, 1/12, 1/12, 1/12,
            3/8, 1/8, 1/8, 3/8,
            3/8, 1/8, 1/12, 1/12, 1/12, 1/4,
            1,
            1/8, 1/8, 1/8, 1/8, 1/8, 1/8, 1/8, 1/8)

writeSong(title = title, key = key, n_pitch = n_pitch, n_durr = n_durr)
```


# Granular Song Conversion

We now use the **recodeSongs()** function to convert our songs into the granular, 1/8 breakdown version. Again, this chunk of code will not be run as all the data was already created and stored in the */recoded_songs/* directory.
```{r convertSongs, eval=FALSE}
songs <- list.files("songs")
for(i in seq_along(songs)){
  df <- read.csv(paste("songs/", songs[i], sep=""), stringsAsFactors = FALSE)
  new_df <- recodeSongs(df)
  print(songs[i])
  print(dim(new_df))
  write.csv(new_df, paste("recoded_songs/", songs[i], sep=""), row.names = FALSE)
}
```


# Model Fitting Approaches

## Approach 1: Full Melody Model

In our first approach we build a model on the full melodic sequence of a song. Each song is considered as an observation fed to our single model, where there is **no assumption** about similarities or differences in melodic structure based on different sections of the song (i.e., Ia, Ib, Ic, Id, IV, V). 

```{r approach1}
# First appending all song melodies to a list
songs <- list.files("songs")
melodies.list <- list()
for(i in seq_along(songs)){
  melodies.list[[i]] <- read.csv(paste("recoded_songs/", songs[i], sep=""), stringsAsFactors = FALSE)$state
}

# Collapsing into one large matrix
melodies.mat <-  matrix(unlist(melodies.list), ncol = length(melodies.list[[1]]), byrow = TRUE)

# Parameters for seqdef() function
alph <- c("C_S", "C_L", "Db_S", "Db_L", "D_S", "D_L",
                       "Eb_S", "Eb_L", "E_S", "E_L", "F_S", "F_L",
                       "Gb_S", "Gb_L", "G_S", "G_L", "Ab_S", "Ab_L",
                       "A_S", "A_L", "Bb_S", "Bb_L", "B_S", "B_L",
                       "R_S", "R_L")
colors <- rep(c(brewer.pal(12, 'Set3'), 'black'), each = 2)

# Creating a TraMineR sequence object used to feed into the pstree() later on.
# Must specify the alphabet (all possible states) as well as cpal which
# is for coloring states in the plotted tree later. 
seq <- seqdef(melodies.mat, alphabet =alph, cpal = colors)

# ~~~~ Fitting model ~~~~~
# For generating
# Change parameters for this to produce different models
pst.g <- pstree(seq, L = 16, nmin=2, ymin=0)

# For classification
# Currently not being used
pst.c <- pstree(seq, L = 16, nmin=2, ymin=0.001)
pst.c <- prune(pst.c, gain="G1", C=1.2)

```


## Approach 2: Model Per Section

In the second approach, we build a single model for each section (i.e., Ia, Ib, IV, Ic, V, and Id) totalling **six models**. We will then assess the similarity of the models through *pairwise divergence* metrics that is analogous to *KL-divergence*. Doing may allow us to identify sections whose melodies repeat, therefore justifying collapsing sections together.
```{r approach2}
songs <- list.files("songs")

# Creating a separate matrix of melodies for each section
Ia <- list()
Ib <- list()
IV <- list()
Ic <- list()
V <- list()
Id <- list()
for(i in seq_along(songs)){
  song <- read.csv(paste("recoded_songs/", songs[i], sep=""), stringsAsFactors = FALSE)
  Ia[[i]] <- song$state[song$section == "Ia"]
  Ib[[i]] <- song$state[song$section == "Ib"]
  IV[[i]] <- song$state[song$section == "IV"]
  Ic[[i]] <- song$state[song$section == "Ic"]
  V[[i]] <- song$state[song$section == "V"]
  Id[[i]] <- song$state[song$section == "Id"]
}

Ia <- matrix(unlist(Ia), ncol = 16, byrow = TRUE)
Ib <- matrix(unlist(Ib), ncol = 16, byrow = TRUE)
IV <- matrix(unlist(IV), ncol = 16, byrow = TRUE)
Ic <- matrix(unlist(Ic), ncol = 16, byrow = TRUE)
V <- matrix(unlist(V), ncol = 16, byrow = TRUE)
Id <- matrix(unlist(Id), ncol = 16, byrow = TRUE)
mel_secs <- rbind(Ia, Ib, IV, Ic, V, Id)
section <- c(rep("Ia", length(songs)),
             rep("Ib", length(songs)),
             rep("IV", length(songs)),
             rep("Ic", length(songs)),
             rep("V", length(songs)),
             rep("Id", length(songs)))

# As a group
seq <- seqdef(mel_secs, alphabet = alph, cpal = colors)

# ~~~~~ Fitting models ~~~~~
# Generating model
# Change parameters for this to produce different models
sections.pst.g <- pstree(seq, L = 8, nmin = 2, ymin = 0.001, group = section)
sections.pst.g <- prune(sections.pst.g, gain = "G1", C = 1.2)

# Classification model
# Currently not being used
sections.pst.c <- pstree(seq, L = 8, nmin = 2, ymin = 0.001, group = section)
sections.pst.c <- prune(sections.pst.c, gain = "G1", C = 1.2)

```

# Pairwise Divergence Assessment
```{r pairwiseDivergence}
# Building models
Ia.seq <- seqdef(Ia, alphabet = alph, cpal = colors)
Ia.pst <- pstree(Ia.seq, L = 8, nmin=2, ymin=0.001)

Ib.seq <- seqdef(Ib, alphabet = alph, cpal = colors)
Ib.pst <- pstree(Ib.seq, L = 8, nmin=2, ymin = 0.001)

IV.seq <- seqdef(IV, alphabet = alph, cpal = colors)
IV.pst <- pstree(IV.seq, L = 8, nmin=2, ymin=0.001)

Ic.seq <- seqdef(Ic, alphabet = alph, cpal = colors)
Ic.pst <- pstree(Ic.seq, L = 8, nmin=2, ymin=0.001)

V.seq <- seqdef(V, alphabet = alph, cpal = colors)
V.pst <- pstree(V.seq, L = 8, nmin=2, ymin=0.001)

Id.seq <- seqdef(Id, alphabet = alph, cpal = colors)
Id.pst <- pstree(Id.seq, L = 8, nmin=2, ymin=0.001)

# Creating dataframe of all combinations and their divergence metric
combos <- combn(1:6, 2) %>% as.matrix() %>% t %>% as.data.frame()
combos$val <- 0
sec.models <- list(Ia.pst, Ib.pst, IV.pst,
                   Ic.pst, V.pst, Id.pst)
for(i in 1:nrow(combos)){
  sub1 <- sec.models[[combos$V1[i]]]

  sub2 <- sec.models[[combos$V2[i]]]

  combos$val[i] <- pdist(sub1, sub2,
                         l = 8, method = "cp",
                         symetric = TRUE, output = "mean",
                         ns = 100)
}

# Recoding sections
combos$V1 <- recode_factor(combos$V1, `1` = "Ia", `2` = "Ib",
                    `3` = "IV", `4` = "Ic", `5` = "V", `6` = "Id",
                    .ordered = TRUE)
combos$V2 <- recode_factor(combos$V2, `1` = "Ia", `2` = "Ib",
                    `3` = "IV", `4` = "Ic", `5` = "V", `6` = "Id",
                    .ordered = TRUE)
# Plotting heatmap
ggplot(data = combos) + 
  geom_tile(mapping=aes(x = V1, y = V2, fill = val)) + 
  geom_text(mapping=aes(x = V1, y = V2, label = round(val,3)),
            colour = "white") + 
  labs(x = "", y = "", title = "Pairwise Divergences of Blues Sections") + 
  scale_fill_continuous(name = "Divergence")

```

# Saving Workspace With Models
```{r saveModels}
keep <- c("convertSeqToNotes", 'pst.g', 'pst.c', 'sections.pst.g', 'sections.pst.c')
rm(list = setdiff(ls(), keep))
save.image("models.RData")
```